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1 .  INTRODUCTION 

The  application  of  finite-difference  methods  in  numerically  solving 
partial  differential  equations  governing  fluid  flow  has  become  increas¬ 
ingly  commonplace  over  the  past  two  decades.  Early  work  was  concentrat¬ 
ed  on  solving  simple  linearized  equations  for  very  simple  geometric  con¬ 
figurations.  As  computers  became  more  sophisticated,  algorithms  were 
improved  and  applications  to  more  complex  problems  were  attempted.  To¬ 
day,  most  companies  and  government  agencies  with  interests  in  fluid  dy¬ 
namics  have  in-house  capabilities  for  solving  flow  problems  using  numer¬ 
ical  techniques. 

Even  though  complicated  problems  can  be  solved,  there  are  still  a 
number  of  pacing  areas  that  are  crucial  in  the  computational  fluid  dy¬ 
namics  field.  One  of  these  areas,  mesh  generation,  is  perhaps  one  of 
the  most  important  topics  needing  further  development  if  continued  prog¬ 
ress  is  to  be  made.  In  fact,  coordinate  system  selection  and  grid  gen¬ 
eration  are  probably  the  most  important  topics  requiring  study  if  con¬ 
tinued  advances  in  digital  simulation  of  fluid  flow  around  flight 
vehicles  are  to  be  made. 

While  a  number  of  problems  in  grid  generation  for  different  geome¬ 
tries  can  be  identified,  the  placement  of  grid  points  in  order  to  ade¬ 
quately  resolve  a  flow  and  provide  a  reduction  of  global  error  in  a  nu¬ 
merical  calculation  is  of  major  interest.  Since  the  solution  of  a 
particular  flow  is  not  known  a  priori,  the  grid  points  cannot  be  placed 
in  the  best  positions  before  the  calculation  is  complete.  Consequently, 
it  is  advantageous  to  adjust  the  grid  point  locations  in  such  a  way  as 
to  provide  the  best  solution  as  the  computation  progresses.  This  idea 
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of  adaptive  grids  was  the  subject  of  AFOSR  Grant-83-0167,  "Application 
of  Adaptive  Grids  in  Solving  the  Partial  Differential  Equations  Govern¬ 
ing  Fluid  Flow." 

The  original  proposal  for  work  under  this  program  was  based  upon  ap¬ 
plying  the  equidistribution  concept  for  generating  grids  to  two-  and 
three-dimensional  problems.  The  technique  applied  by  Rai  and  Ander¬ 
son1'2  was  selected  as  a  candidate  method.  This  method  was  used  in  ear- 
.y  experiments  in  two-dimensions  for  small  amounts  of  grid  adjustment 
with  good  success. 

Early  in  the  program,  several  example  calculations  were  completed  for 
high  values  of  grid  adaption  with  the  Rai  and  Anderson  scheme.  For 
large  adaption  rates,  severe  grid  skewness  was  encountered.  The  one-di¬ 
mensional  method  employed  by  Dwyer  et  al.3  was  also  extended  to  two  di¬ 
mensions.  Similar  examples  were  tested  with  this  scheme  and  severe  grid 
skewness  also  resulted  for  large  adaption  rates.  Faced  with  severe  grid 
distortion  in  the  multidimensional  case,  a  new  way  of  creating  an  adap¬ 
tive  mesh  was  needed.  Not  only  must  the  technique  be  capable  of  provid¬ 
ing  the  desired  grid  adaption  but  some  positive  method  of  grid  skewness 
control  is  necessary.  It  is  now  accepted  by  researchers  that  any  method 
for  generating  an  adaptive  grid  without  an  active  skewness  control  will 
ultimately  fail  in  two  or  three  dimensions. 

Faced  with  this  mesh  control  problem,  it  was  necessary  to  review  the 
fundamental  concept  of  equidistribution  of  a  weight  function  over  an 
area.  This  concept  was  retained  as  the  main  idea  in  the  generation  of 
an  adaptive  grid.  However,  it  was  clear  that  the  relationship  between 
the  weight  function  and  the  mesh  skewness  must  be  developed.  Research 
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on  this  issue  progressed  until  an  explicit  formulation  relating  the  an¬ 
gle  between  mesh  lines  to  the  weight  function  was  developed.  This  also 
lead  to  a  simple  way  of  generating  an  orthogonal  grid  for  two  dimen¬ 
sions  . 
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2.  SIGNIFICANT  MILESTONES  ACHIEVED  AND  STATUS  OF  RESEARCH 
UNDER  AFOSR  GRANT-83-0167 


2.1  Milestones 

During  the  present  program,  a  number  of  important  contributions  were 
made  in  the  understanding  of  adaptive  grids. 

1.  In  an  invited  review  paper  presented  to  ASME , *  a  number  of 
methods  for  creating  two-dimensional  mesh  systems  were  shown  to  be 
virtually  the  same.  These  schemes  are  based  upon  the  concept  of 
equidistribution  of  a  weight  function  over  an  area. 

2.  At  the  same  time  it  was  shown  that  all  schemes  for  generating 
an  adaptive  mesh  without  an  active  skewness  control  will  fail. 

3.  Recently,  the  direct  relationship  between  the  skewness  of  mesh 
lines  and  the  weight  function  was  demonstrated. 5  This  shows  how 
the  grid  distortion  can  be  controlled  while  still  providing  an 
adaptive  grid. 

2.2  Status 

At  the  time  of  the  termination  of  this  program,  the  concept  of  con¬ 
structing  a  two-dimensional  grid  based  on  equidistribution  was  well  es¬ 
tablished.  In  addition,  the  coupling  between  grid  adaption,  skewness, 
and  weight  function  was  well  understood.  A  series  of  numerical  experi¬ 
ments  is  needed  to  demonstrate  that  grid  adaption  and  control  can  both 
be  attained  for  two-dimensional  grids.  This  demonstration  would  pave 
the  way  for  practical  application  to  the  equations  governing  fluid  flow 
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In  reviewing  the  goals  and  proposed  methodology  for  attaining  the 
goals  of  this  program,  few  changes  were  required.  The  necessary  changes 
in  approach  were  a  result  of  knowledge  gained  as  research  results  were 
obtained.  Continued  research  on  adaptive  grids  using  the  ideas  devel¬ 
oped  on  this  program  should  provide  a  practical  method  that  can  be  ap¬ 


plied  to  useful  problems. 
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Fig.  9  Grid  adaption  for  sinusoidal  function 
distribution,  a  -  0,  b  -  10 
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Fig.  11  Comparison  of  error  in  surface  potential 
for  a  cylinder 


Fig.  10  Grid  adaption  for  sinusoidal  function 
distribution,  a  -  2,  b  »  3 


Fig.  12  Grid  structure  for  Mach  reflection  of 
shock  from  a  ramp 
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Fig.  4  Early  time  grid  structure  and  isotherms 
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of  grid  smoothness  is  invoked,  an  approximate  expression 
for  the  error  terms  can  be  derived  in  terms  of  the  co¬ 
efficient  matrices  of  the  original  system  and  the  errors 
incurred  in  forming  the  first  derivatives  of  the  depen¬ 
dent  variables.  These  error  expressions  are  weighted 
averages  of  all  the  contributions  of  the  governing 
equations  for  the  fluid  flow  problem. 

AREAS  OF  CONTINUED  AND  FUTURE  RESEARCH 

A  number  of  problem  areas  need  considerable  work 
in  order  to  create  adaptive  grid  schemes  which  jure 
readily  applicable  to  a  wide  class  of  problems.  These 
areas  include  the  following: 

1.  Better  estimates  of  numerical  error  are  nec¬ 
essary  for  use  in  those  adaptive  methods  which  attempt 
to  reduce  the  error  in  the  solution. 

2.  Additional  effort  must  be  expended  in  gaining 
an  understanding  of  how  points  are  allocated  in  dif¬ 
ferent  regions  for  two-  and  three-dimensional  problems 
for  a  given  grid  adaption  scheme. 

3.  New  methods  for  controlling  point  motion  while 
still  providing  sufficient  adaption  should  be  developed 
for  the  various  schemes. 

4.  Application  of  existing  adaptive  methods  to 
three-dimensional  problems  should  be  encouraged  because 
the  largest  gains  in  computational  efficiency  will  be 
realized  in  these  cases. 

Research  in  adaptive  grid  schemes  has  really 
received  emphasis  during  approximately  the  past  five 
years.  Fresh  and  innovative  ideas  are  needed  for  con¬ 
struction  of  better  approaches  to  the  problems  in  this 
area.  Researchers  should  be  encouraged  to  thoroughly 
explore  any  idea  that  may  have  an  impact  on  the  devel¬ 
opment  of  new,  successful  techniques. 

ACKNOWIEDGEMENTS 

This  work  was  partially  sponsored  by  funds  made 
available  through  U.S.  Army  Contract  DAAK11-82-R-0123 
and  NASA  Grant  NC,T  016-002-801. 

Figures  4,  5,  and  14  through  20  from  cited  refer¬ 
ences  are  reprinted  with  permission  of  the  American 
Institute  of  Aeronautics  and  Astronautics. 

Figures  2  and  3  are  reprinted  by  permission  of  the 
publisher  from  the  cited  article  in  Numerical  Grid 
Generation ,  Joe  F.  Thompson  (ed.).  North  Holland  1982, 
copyright  1982  by  Elsevier  Science  Publishing  Co.,  Inc. 


REFERENCES 

1.  Ablow,  C.M.,  "Equidistant  Mesh  for  Gas  Dynamic 
Calculations,"  Numerical  Grid  Generation,  Joe  F.  Thomp¬ 
son,  editor,  Elsevier  Science  Publishing  Co.,  Inc., 
1982,  pp.  859-863. 

2.  Acharya,  S.  and  Patankar,  S.V.,  "Use  of  An 
Adaptive  Grid  for  Parabolic  Flows,"  AIAA  Paper  82-1015, 
presented  at  the  AIAA/ASME  3rd  Joint  Thermophysics, 
Fluids,  Plasma  and  Heat  Transfer  Conference,  St.  Louis, 
Missouri,  June  1982. 


3.  Anderson,  D.A, ,  “Solution  Adaptive  Grids  for 
Partial  Differential  Equations,"  ARO  Report  82-3, 
Proceedings  of  the  1982  Army  Numerical  Analysis  and 
Computers  Conference,  Vicksburg,  Mississippi,  February 
1982,  pp.  575-591. 

4.  Anderson,  D.A.  and  Rai,  M.M. ,  "A  New  Approach 

to  Solution  Adaptive  Grids,"  Proceedings  of  the  AIAA/ASME 
Symposium  on  Computers  in  Flow  Predictions  and  Fluid 
Dynamics  Experiments,  November  1981,  p.  95. 

5.  Anderson,  D.A.  and  Rai,  M.M. ,  "The  Use  of 
Solution  Adaptive  Grids  in  Solving  Partial  Differential 
Equations,"  Numerical  Grid  Generation,  Joe  F.  Thompson, 
editor,  Elsevier  Science  Publishing  Co.,  Inc.,  1982, 
pp.  317-338. 

6.  Brackbill,  J.U. ,  "Coordinate  System  Control: 
Adaptive  Meshes,"  Numerical  Grid  Generation,  Joe  F. 
Thompson,  editor,  Elsevier  Science  Publishing  Co.,  Inc., 
1982,  pp.  277-294. 

7.  Brackbill,  J.  and  Saltzman,  J. ,  "Adaptive  Zoning 
for  Singular  Problems  in  Two  Dimensions,”  Journal  of 
Computational  Physics,  Vol.  46,  June  1982,  pp.  342-368. 

8.  Dwyer,  H.A.,  "Grid  Adaption  for  Problems  with 
Separation,  Cell  Reynolds  Number,  Shock-Boundary  Layer 
Interaction,  and  Accuracy,"  AIAA  Paper  83-0449,  presented 
at  AIAA  21st  Aerospace  Sciences  Meeting,  Reno,  Nevada, 
January  1983. 

9.  Dwyer,  H.,  Kee,  R. ,  and  Sanders,  B. ,  "Adaptive 
Grid  Method  for  Problems  in  Fluid  Mechanics  and  Heat 
Transfer,"  AIAA  Journal ,  Vol.  18,  October  1980,  pp.  1205- 
1212. 

10.  Dwyer,  H.A. ,  Smooke,  M.D.,  and  Kee,  R.J., 
"Adaptive  Gridding  for  Finite  Difference  Solutions  to 
Heat  and  Mass  Transfer  Problems,"  Numerical  Grid  Geher- 
ation,  Joe  F.  Thompson,  editor,  Elsevier  Science 
Publishing  Co.,  Inc.,  1982,  pp.  339-356. 

11.  Gelinas,  R.J.  and  Doss,  S.K. ,  "The  Moving 
Finite  Element  Method:  Applications  to  General  Partial 
Differential  Equations  with  Multiple  Large  Gradients," 
Journal  of  Computational  Physics,  Vol.  40,  1981, 

pp.  202-249. 

12.  Gnoffo,  P.A. ,  "A  Vectorized,  Finite-Volume, 
Adaptive-Grid  Algorithm  for  Navier-Stokes  Calculations," 
Numerical  Grid  Generation,  Joe  F.  Thompson,  editor, 
Elsevier  Science  Publishing  Co.,  Inc.,  1982,  pp.  819-835. 

13.  Gnoffo,  P.A. ,  "A  Vectorized,  Finite-Volume, 
Adaptive  Grid  Algorithm  Applied  to  Planetary  Entry 
Problems,"  AIAA  Paper  82-1018,  presented  at  AIAA/ASME 
3rd  Joint  Thermophysics,  Fluids,  Plasma  and  Heat  Trans¬ 
fer  Conference,  St.  Louis,  Missouri,  June  1982. 

14.  Hindman,  R.G.,  Kutler,  P.,  and  Anderson,  D.  , 
"TVo-Dimensional  Unsteady  Euler-Equation  Solver  for 
Arbitrarily  Shaped  Flow  Regions,”  AIAA  Journal,  Vol.  19, 
April  1981,  pp.  424-431. 

15.  Hindman,  R.G.  and  Spencer,  J. ,  "A  New  Approach 
to  Truly  Adaptive  Grid  Generation,"  AIAA  Paper  83-0450, 
presented  at  AIAA  21st  Aerospace  Sciences  Meeting,  Reno, 
Nevada,  January  1983. 


GRID  POINT  CONTROL  AND  GRID  ADAPTION 


0n«»  of  the  most  difficult  problems  encountered  in 
using  adaptive  grids  is  that  of  grid  point  control. 

This  problem  is  not  severe  in  one-dimensional  adaption 
but  is  much  more  acute  in  multidimensional  applications. 
In  the  variational  formulation  of  Ref.  7,  grid  control 
was  a  consideration  in  the  construction  of  the  method. 
The  problem  remains  in  selecting  the  values  of  X  that 
are  used  to  provide  the  contribut ions  due  to  smoothness, 
orthogonality,  and  adaption.  In  order  to  provide  suf¬ 
ficient  adaption,  high  grid  skewness  may  be  encountered. 
Computed  results  show  that  a  direct  trade-off  exists 
among  the  various  terms  in  the  grid  generator.  Once 
again,  the  problem  of  selecting  certain  parameters 
before  the  solution  is  obtained  occurs. 

In  Dwyer's  (8)  recent  paper,  he  describes  a  tech¬ 
nique  of  predetermining  the  percentage  of  mesh  points 
assigned  to  grid  adaption.  In  this  paper,  a  one¬ 
dimensional  adaption  technique  was  used  similar  to  that 
in  Eq.  (29).  The  weight  function  used  was  of  the  form 


what  is  the  best  choice  of  the  weighting  function  in  the 
variational  schemes  and  what  is  the  best  choice  for  e 
in  Eq.  (57)?  This  question  is  not  easily  answered  even 
when  a  single  scalar  equation  is  used. 

In  viscous  problems,  resolution  of  viscous  regions 
is  probably  best  accomplished  by  keeping  the  cell 
Reynolds  number  less  than  one.  This  provides  a  well 
defined  viscous  layer  where  first-order  upwind  schemes 
will  not  produce  large  artificial  diffusion  and  an 
oscillation  free  solution  will  be  obtained  when  higher- 
order  schemes  are  used.  In  principle,  any  function 
which  would  permit  clustering,  such  as  gradient  infor¬ 
mation  should  provide  adequate  adaption.  Numerous 
researchers  have  used  second  derivative  information. 

In  Ref.  8  a  nonlinear  combination  of  first  and  second 
derivatives  has  been  used  while  second  derivatives  have 
also  been  used  in  Refs.  5  and  21.  White  (28)  has  used 
curvature  of  the  solution  for  the  clustering  function 
and  has  shown  good  results  for  a  scalar  equation.  It 
should  be  noted  that  a  second  derivative  clustering 
function  can  be  viewed  as  a  first  approximation  to  the 
curvature . 


=  1  +  b 


|5li 


(60) 


where  b  is  a  constant  and  f  is  the  function  which  is 
monitored  and  used  for  adaption.  The  ratio,  ,  is 
formed  by  computing  the  relative  contribution  of  the 
grid  adaption  to  the  computational  coordinate  and  is 
given  by 
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If  the  ratio  is  held  at  a  fixed  value  during  a  cal¬ 
culation  and  the  value  of  b  is  determined  from  Eq.  (61) , 
the  relative  weight  placed  on  adaptivity  in  the  mesh 
remains  fixed.  This  mechanism  provides  some  control  on 
the  assignment  of  mesh  points  to  different  areas  of 
interest  in  the  problem.  A  similar  division  can  be 
used  in  conjunction  with  the  inverse  adaption  scheme 
Of  Eq.  (24). 


Grid  point  control  for  other  schemes  such  as  that 
of  Rai  and  Anderson  is  achieved  by  logic  built  into  the 
adaption  algorithm  as  opposed  to  the  governing  equations. 
For  example  the  constant  K  in  the  one-dimensional  grid 
speed  equation  (Eq.  (57)  |  is  continuously  adjusted 
during  the  numerical  computations  to  prevent  the  maximum 
grid  speed  from  exceeding  a  predetermined  value.  An 
additional  constraint  is  placed  upon  the  grid  speeds 
when  the  grid  points  are  closer  than  a  specified  value. 

To  prevent  excessive  stretching,  the  grid  speeds  at 
those  points  are  attenuated  by  the  factor 


-i/x? 


cal 


where  is  the  grid  speed  computed  from  the  grid 

cai 

speed  equation.  This  prevents  grid  speeds  from  being 
excessively  large  in  regions  where  point  density  is  high 
and  provides  good  control  of  the  grid  point  motion. 


Another  question  closely  associated  with  grid  point 
control  is  that  of  attempting  to  define  appropriate 
functions  to  use  in  the  adaption  process.  For  example, 


For  each  grid  adaption  scheme  and  each  problem,  the 
choice  of  a  grid  adaption  function  may  be  different. 

For  example,  the  choice  of  gradient  alone  for  the  weight 
function  w^  in  Eq.  (22)  would  not  be  satisfactory . 
Methods  based  upon  this  expression  would  produce  grids 
with  no  mesh  points  in  uniform  regions.  On  the  other 
hand,  gradient  may  work  well  in  the  two-dimensional 
application  of  Eq.  (18)  because  measures  of  grid  smooth¬ 
ness  and  orthogonality  are  also  included. 

As  previously  noted,  the  lowest-order  error  term 
of  the  modified  partial  differential  equation  was  used 
to  drive  the  grid  in  Refs.  20,  21,  and  22.  This  error 
term  involves  a  third  derivative  assuming  a  second-order 
method  is  used.  Numerical  evaluation  of  the  third 
derivative  (or  any  higher-order  derivative)  is  generally 
very  noisy.  As  a  result,  data  used  to  evaluate  deriv¬ 
atives  is  smoothed  to  prevent  feedback  into  the  grid 
motion.  If  the  assumption  of  a  very  smooth  grid  is 
made,  the  error  can  frequently  be  approximated  with 
lower-order  derivatives. 


When  systems  of  equations  are  solved  or  physical 
problems  are  studied  where  the  resolution  of  more  than 
one  physical  event  is  necessary,  additional  difficulties 
are  encountered.  In  Ref.  8,  the  problem  of  flame  propa¬ 
gation  about  a  spherical  particle  in  the  presence  of  a 
low  Reynolds  number  Stokes  flow  was  studied.  In  this 
case,  resolution  of  both  the  viscous  regions  and  the 
laminar  flame  zone  is  desired.  The  grid  adaption  func¬ 
tions  should  include  contributions  which  would  resolve 
both  the  flame  front  and  the  viscous  region.  In  this 
case,  the  problem  was  solved  on  two  grids  and  the  solu¬ 
tion  was  determined  by  interpolating  between  them.  The 
two  grids  showing  details  of  each  region  are  shown  in 
Figs.  19  and  20. 

A  one-dimensional  inviscid  shock  tube  problem  was 
studied  in  Ref.  16.  in  this  example,  a  reduction  of 
the  truncation  error  in  the  modified  differential  equa¬ 
tion  was  desired.  An  adaptive  grid  was  used  and  the 
grid  driving  function  (wJ)  was  selected  to  be  a  linear 
combination  of  the  error  terms  from  each  of  the  governing 
equations.  A  more  analytical  approach  to  the  deriva¬ 
tion  of  the  function  used  to  drive  the  adaptive  grid 
was  taken  in  Ref.  5.  The  Euler  equations  in  multidimen¬ 
sional  space  were  under  study  and  expressions  for  the 
error  terms  in  Eq.  (58)  were  needed.  When  the  assumption 


a  rectangle  in  computational  space.  Since  the  grid  is 
three  sided  in  physical  space,  a  geometric  singularity 
is  shown  on  the  shock  wave.  It  should  be  remembered 
that  the  grid  point  locations  for  this  scheme  are  com¬ 
puted  by  integrating  the  grid  speeds  instead  of  solving 
the  steady  grid  equation. 

Hindman  and  Spencer  (15)  have  continued  this  ap¬ 
proach  and  have  considered  the  one-dimensional  grid 
equat ion 

3_ 

Xf,C  +  X£P  “  0  <54) 

Again,  the  grid  speeds  are  established  by  differentiating 
this  expression  with  respect  to  T.  The  P  function  was 
selected  (at  least  for  one  case)  so  that  Eq .  (51)  was 
consistent  with  Eq.  (42).  The  form  for  P  becomes 


Figure  14  shows  the  time  history  of  the  grid  motion  when 
the  inviscid  Burgers'  equation  is  solved  using  a  single 
discontinuity  for  initial  conditions.  MacCormack 'a 
scheme  was  used  to  integrate  both  the  equations  of 
motion  and  the  grid  speeds.  The  main  result  of  interest 
in  this  figure  is  that  the  grid  tends  to  relax  as  time 
increases.  This  can  be  avoided  by  solving  the  steady 
grid  equation  after  a  predetermined  number  of  steps  and 
interpolating  the  solution  on  the  new  grid. 

Rai  and  Anderson  (20,21,22)  and  Anderson  and  Rai 
(3,4,5)  have  constructed  an  adaptive  grid  scheme  using 
a  different  approach  from  those  discussed  above.  This 
method  also  determines  the  grid  speed  and  the  grid 
point  locations  are  established  by  integration.  A 
much  simpler  grid  speed  expression  than  that  used  in 
Ref.  14  was  developed.  The  method  is  based  upon  an 
attraction  model.  It  is  assumed  that  the  bast  grid  for 
a  given  problem  is  one  where  the  solution  error  at 
every  point  reaches  a  constant  value.  In  creating  this 
grid,  more  points  are  needed  in  regions  where  the  error 
is  larger  than  the  average  and  fewer  points  are  needed 
in  regions  where  the  local  error  is  smaller  than  the 
average  value.  Tins  leads  directly  to  the  idea  of 
associating  a  capacity  to  induce  velocity  with  the  local 
error  at  each  point.  For  two  points,  A  and  B,  the  grid 
speed  induced  at  point  B  due  to  an  error  at  point  A  is 
written 


V  =  K 
BA 


UIA  -  Hav 


whore  c  represents  the  local  error,  av  indicates  average 
value  over  all  points,  raA  is  the  distance  from  A  to  B 
in  computational  space,  K  is  a  proport ional ity  constant, 
and  n  is  a  power  which  controls  the  attenuation  of  the 
attraction  with  distance.  For  a  one-dimensional  problem, 
the  grid  speed  in  physical  space  becomes 
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Grid  point  positions  produced  by  integrating  Eq.  (57) 
must  not  cross.  In  this  scheme,  grid  points  do  not 
cross  for  sufficiently  small  time  steps  because  the  local 
value  of  the  error  will  be  less  than  the  average  value 
if  points  are  very  close  together.  This  creates  a  sign 
switch  on  the  grid  speed  (repulsion).  As  points  approach 
each  other  the  term  £x  becomes  very  large  which  makes 
further  movement  of  the  grid  points  very  difficult.  The 
constant,  K,  in  Eq.  (57)  is  adjusted  to  scale  the  grid 
speed  to  a  predetermined  maximum  value  during  the  calcu¬ 
lations  . 

7Tie  driving  function  used  to  establish  grid  point 
motion  using  this  scheme  (|e|  -  |e|  )  can  be  based  on 

error  or  the  gradient  of  any  flow  variable  or  any  other 
function  which  provides  the  desired  attraction.  A  num¬ 
ber  of  error  measures  were  used  to  drive  the  grids  in 
Refs.  5  and  21.  In  general,  the  form  of  the  error  was 
established  by  using  an  approximation  to  the  lowest- 
order  error  term  of  the  modified  partial  differential 
equation.  Thompson  (24)  suggests  that  the  grid  solution 
obtained  using  Eq.  (57)  is  equivalent  to  solving  a 
variational  problem  by  an  iterative  approach,  in  this 


case,  the  variation  of 


is  minimized  over  the 


field.  This  corresponds  to  the  weight  function  w  in 
Eq.  (22).  1 

In  order  to  provide  a  smooth  grid  at  the  boundaries, 
a  reflection  at  the  boundaries  in  the  computational 
domain  was  used.  For  example,  at  the  right-hand  bound¬ 
ary,  £  “  1  and  £  =  0.  However,  if  this  condition  is 

explicitly  used,  sometimes  the  grid  will  not  be  smooth 
near  this  boundary.  If  a  series  of  points  are  reflected 
and  have  an  error  measure  assigned  to  them  in  such  a  way 
that  /£.\  “  0  when  computed  from  the  grid  equation, 

v  /e=i 

the  boundary  region  will  have  a  very  smooth  grid.  The 
Influence  of  this  reflection  is  in  the  value  of  along 
the  boundary.  This  value  and  the  resulting  point  loca¬ 
tions  are  very  sensitive  to  boundary  point  treatment. 

Figure  15  shows  the  converged  adaptive  grid  gener¬ 
ated  for  the  supersonic  blunt  body  problem  in  two 
dimensions.  Figure  16  shows  the  difference  between  the 
total  enthalpy  at  the  cylinder  surface  as  computed  from 
the  numerical  solution  and  the  free  stream  value.  Since 
this  is  an  inviscid  calculation,  this  is  a  measure  of 
the  solution  accuracy.  Error  in  total  enthalpy  is 
significantly  reduced  when  an  adaptive  grid  is  used. 

A  technique  for  aligning  a  grid  with  a  high  gradi¬ 
ent  region  is  presented  in  Refs.  5  and  22.  This  tech¬ 
nique  was  developed  for  use  with  shock  capturing  methods. 
It  is  well  known  that  shock  aligned  coordinate  systems 
permit  much  better  computation  of  shocks  with  these 
methods  because  the  flux  terms  are  then  continuous 
across  the  shocks,  shock  alignment  is  accomplished  by 
creating  a  grid  speed  in  a  two-dimensional  problem  along 
only  one  of  the  coordinates.  This  grid  speed  is  propor¬ 
tional  to  the  product  of  the  gradient  of  some  property 
of  the  solution  (density,  pressure)  along  both  the  f,  and 
H  directions.  The  result  is  an  effective  rotation  of 
the  coordinate  line  segments  in  such  a  way  as  to  line  up 
with  level  surfaces  of  h.  Again,  only  movement  in  one 
coordinate  direction  is  necessary  in  a  two-dimensional 
problem.  Figure  17  shows  the  converged  grid  for  an 
oblique  shock  and  the  pressure  distribution  for  both 
shock  capturing  on  a  uniform  grid  and  an  aligning  grid 
is  shown  in  Fig.  18.  The  quality  of  the  solution  is 
dramatically  improved  when  the  aligning  grid  is  used. 


(57) 


where  it  is  assumed  that 


Iq^lAC  -  2g  >  0 

In  regions  where  this  is  not  true,  any  positive  value 
of  a  satisfies  the  condition  on  grid  crossing. 

In  Ref.  17  an  adaptive  grid  generation  scheme  is 
proposed  which  for  one-dimensional  problems  is  of  the 
form 


2 

xf;f  +  ATjr*£  “  0  (44) 

where  A  is  a  positive  constant  and  T  is  some  nonnegative 
function  of  the  solution.  In  the  context  of  Eq.  (22), 
this  can  be  interpreted  as  a  scheme  which  is  similar  to 
either  Dwyer's  or  Gnoffo's  method  with  the  weight  func¬ 
tion  defined  as 


E 

/  AT^x^df, 


(45) 


The  two-dimensional  formulation  of  Eq.  (44)  is  of  the 
form 


All  of  the  methods  discussed  thus  far  have  relied 
upon  solving  a  steady  grid  equation  to  determine  the 
grid  point  location.  The  grid  speed  is  determined  by 
using  a  backward  difference  using  the  computed  mesh 
point  location  and  the  previously  known  positions. 
Hindman  et  al.  (14)  developed  a  technique  of  computing 
the  grid  speed  directly  from  the  grid  generator.  The 
idea  is  to  evaluate  the  time  derivative  of  the  steady 
grid  generation  equation  and  solve  this  equation  for  the 
grid  speeds.  In  principle,  this  idea  is  applicable  to 
the  methods  for  creating  an  adaptive  grid  presented 
above. 


If  the  Thompson  scheme  is  used  as  a  starting  point, 
the  equations  which  determine  the  mapping  relating 
physical  and  computational  space  are 


V2?  -  P 

T2  n  -  q 


(49) 


In  the  computational  domain  these  equations  became 


G[x]  -  0 
G  [y  ]  ■  0 


(50) 
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The  equations  which  are  solved  in  the  computational 
domain  are 


^CE  '  26xen  +  Yxnn  +  J  [asExE  +  BVn'  *  0 


^EE  -  2ByEn  +  YYnn  +  J  [ASEyE  +  BVn1  *  0 


(47) 


where  the  coefficients  a,  8,  and  y  are  functions  of  the 
metrics,  and  S  and  T  are  the  grid  adaption  functions. 

As  in  the  one-dimensional  case,  an  estimate  of  the  up¬ 
per  limit  for  A  and  B  can  be  established  by  requiring 
that  grid  lines  be  noncrossing.  These  estimates  are 


a  <  2/(y|sf|Af,) 
b  <  2/(a|T  |dn) 


(48) 


The  potential  flow  ak>out  a  cylinder  was  computed 
in  Ref.  17  as  an  example  using  the  adaptive  grid  scheme 
in  Kq.  (4f>)  .  Figure  11  shows  the  error  between  the 
exact  and  computed  solutions  for  the  surface  potential 
for  a  nonadaptive  conformal  grid  and  an  adaptive  grid* 

In  general,  an  adaptive  grid  provides  solutions  which 
exhibit  lower  error.  In  this  case  the  S  and  T  functions 
were  defined  as 


S  *  (<t>  -  x)^f 


T 
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nn 


and  $  is  the  disturbance  potential. 


where  the  operator  G  is  defined  by 


G 


2fP  — 
p  n 
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and  a,  6,  and  Y  are  the  usual  functions  of  the  metrics. 
Solution  of  Eq.  (46)  provides  the  grid  point  locations 
at  any  time.  The  grid  point  speeds  of  boundary  points 
are  usually  obtained  from  shock  relations  or  other 
expressions  which  must  be  satisfied.  The  interior  grid 
speeds  (xT,  yT)  could  be  obtained  from  backward  differ¬ 
ences  as  noted  previously,  but  a  better  approach  may  be 
to  differentiate  Eqs.  (50)  with  respect  to  T.  This 
yields  a  system  of  equations  of  the  form 

[M]2  -  r  (52) 


where 


z  *  <v  vT 

and 

r  -  -J2(Ptx^  +  QTxn,  PTy^  +  0Tyn>T  (53) 

The  solution  of  Eq.  (52)  yields  the  interior  point  grid 
speed  necessary  to  advance  the  governing  equations  of 
the  fluid  flow  problem.  If  P  and  g  are  selected  tc  be 
functions  of  the  solution  of. the  flow  equations,  the 
grid  adjusts  adaptively  through  the  time  derivatives  of 
these  terms.  The  boundary  grid  point  speeds  influence 
the  interior  solution  through  the  boundary  conditions 
on  the  system  of  equations. 

In  Ref.  14,  the  P  and  g  functions  were  set  equal 
to  rero  providing  grid  and  grid  speed  solutions  which 
did  not  adapt  to  internal  flow  changes  but  only  to  the 
boundary  motion.  Figure  12  shows  the  geometry  of  the 
grid  produced  for  solving  the  problem  of  a  planar  shock 
wave  passing  over  an  inclined  ramp.  Figure  13  shows  the 
grid  used  to  solve  for  the  inviscid  supersonic  flow  over 
an  ogive.  This  grid  in  physical  space  is  mapped  into 
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Replacing  the  derivative  with  respect  to  s  in  Eg.  (31) 
results  in  the  expression 
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This  expression  can  be  expanded  and  becomes  a  second- 
order  differential  equation  in  x  and  y.  The  companion 
relationship  along  constant  f;  surfaces  can  be  derived 
in  the  same  manner. 


may 


The  inverse  relationship 
be  written 
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corresponding  to  Eq.  (29) 


(34) 


The  differential  equation  satisfied  by  S  corresponding 
to  Eq.  (30)  is 


H  (wiV  =  0 


(35) 
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In  this  case,  the  partial  differential  equation 
which  must  be  satisfied  is  easily  recovered  using  the 
fact  that 

sZ*\[*fyl  <3? 

The  governing  expression  becomes 


and  the  grid  adaption  was  implemented  for  a  specified 
function  u(x,y)  on  a  grid  which  is  21  x  21  in  dimension. 

Figure  6  shows  the  results  for  specifying  u(x,y) 
to  be  a  function  of  x  only.  As  expected,  the  grid 
adjusts  only  in  the  x  direction  and  no  changes  occur 
along  the  y  coordinate.  The  clustering  is  achieved 
using  u(x,y)  *  1  for  x  ^  9  and  u(x,y)  -  0  for  x  >_  11. 

The  function,  u(x,y>,  wa3  selected  to  decrease  linearly 
from  one  to  zero  from  x  *  9  to  x  *  11.  .The  two  bands 
of  clustering  which  appear  in  Fig.  6  are  due  to  aver¬ 
aging  used  on  the  grid  point  calculation.  The  aver¬ 
aging  gives  the  effect  of  using  a  second  derivative  for 
clustering  in  this  case.  If  no  averaging  is  used,  the 
grid  clusters  between  x  *  9  and  x  *  11  where  the  deriv¬ 
ative  3u/9x  attains  a  maximum. 

Figure  7  shows  the  grid  resulting  from  applying 
the  adaption  routine  for  a  u(x,y)  which  varies  linearly 
from  one  to  zero  but  this  variation  occurs  about  the 
syauetric  straight  line 

y  -  2.5x  -  15 

passing  through  the  domain.  The  value  of  u  is  taken  to 
be  either  zero  or  one  outside  this  small  two-unit  wide 
region.  In  Fig.  7,  b  »  0  and  the  attraction  is  only 
along  the  n  «  constant  surfaces.  In  this  case,  we  see 
that  the  adaption  in  the  n  direction  also  creates  a 
contraction  along  the  £  -  constant  lines.  This  is  a 
result  of  the  coupling  between  the  two  directions.  A 
case  which  includes  adaption  in  both  directions  is 
shown  in  Fig.  8.  Additional  cases  with  the  correspon¬ 
ding  u(x,y)  are  shown  in  Figs.  9  and  10.  The  function, 
u(x,y) ,  was  set  equal  to  one  above  and  zero  below  a 
sinusoidal  surface  which  shows  up  very  well  in  these 
grid  plots.  This  function  again  decreased  from  one  to 
Zero  over  a  two  Ay  interval  in  the  y  direction.  The 
grid  distortion  problems  inherent  with  this  scheme  are 
apparent  in  these  results.  Notice  that  the  grid  becomes 
distorted  all  along  the  dividing  surface  of  u(x,y)  but 
is  much  worse  where  the  sine  wave  slope  is  roughly  at 
equal  angles  to  the  x  and  y  axes.  At  points  in  this 
neighborhood,  clustering  along  both  directions  creates 
severe  strain  of  the  area  elements  and  leads  to  the 
distortion  shown  in  Figs.  9  and  10. 
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While  the  author  has  not  carried  out  the  algebraic 
manipulations  (lazy)  required  in  Eq.  (33)  ,  it  is  not 
apparent  that  the  result  would  be  equivalent  to  Eq.  (38). 

It  should  be  noted  that  Eq.  (38)  is  the  expression 
which  is  valid  along  an  n  =  constant  surface.  The 
companion  expression  valid  along  a  constant  £  surface 
may  be  written  directly  as 


The  adaption  process  works  well  for  these  cases, 
but  it  is  difficult  to  estimate  the  proper  values  of 
a  and  b  unless  some  numerical  experiments  are  performed. 
It  should  be  noted  that  the  definition  of  the  transfor¬ 
mation  implies  that  arbitrarily  large  values  of  a  and  b 
can  be  used.  Due  to  the  discretization  of  the  differ¬ 
ential  equation,  the  clustering  constants  a  and  b  can 
not  be  arbitrarily  large.  Mastin  and  Thompson  (17)  have 
shown  that  for  a  one-dimensional  problem 

1*^1  <  2|xtl/A?  (41) 
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It  is  of  interest  to  consider  the  grid  adaption 
scheme  resulting  from  the  simultaneous  solution  of  Eqs. 
(38)  and  (39)  and  some  examples  were  studied.  The 
weight  functions  w^  and  v2  were  selected  to  be  of  the 
form 


in  order  to  prevent  the  grid  lines  from  crossing.  For 
the  differential  equation 
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this  leads  to  a  bound  for  the  constant,  a,  which  may  be 
written 


w1  -  1  +  a|u  | 
w2  -  1  +  bluj 
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(40b) 


Many  applications  require  grid  adaption  in  only 
one  dimension.  For  this  reason,  it  is  worthwhile  to 
consider  minimizing  the  functional,  1^,  defined  in 
Eq.  (14)  specialized  to  one  dimension. 
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The  Euler -Lagrange  equation  may  be  written 
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Equation  (24)  is  the  integral  form  of  the  discrete 
adaptive  grid  scheme  proposed  by  Gnoffo  (12,13).  dioffo 
has  used  this  adaptive  grid  generator  while  solving  the 
Navier-Stokes  equations  with  a  weight  function  based 
upon  monitoring  either  Mach  number  or  velocity  gradient, 
i.e. 


w1  «  1  +  ag 

where 
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In  this  case,  a  first  integral  may  be  directly  written 
as 
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(22) 


This  expression  states  that  the  product  of  the  mesh 
spacing  and  the  weight  function,  w^,  should  remain  con¬ 
stant  in  physical  space.  Equation  (22)  may  be  inte¬ 
grated  to  obtain  the  expression  for  either  the  physical 
coordinate  in  terms  of  the  computational  coordinate  or 
vice  versa.  Let  x  -  0  when  £  *  0  and  x  *  Xgiax  w^en 
£  «  f,max-  If  Eq.  (22)  is  integrated  and  the  computa¬ 
tional  coordinate  is  determined. 
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and  if  the  physical  coordinate  is  evaluated  one  obtains 
the  equation 


x  »  x 


max 


(24) 


Equation  (23)  is  exactly  the  expression  used  by 
Dwyer  (8)  and  Dwyer  et  al.  (9,10).  Many  of  the  appli¬ 
cations  of  this  law  for  grid  adjustment  have  been  in 
combustion  and  heat  or  mass  transfer  and  the  results 
have  been  very  good.  In  Dwyer's  formulation,  the 
weighting  function  w^  was  selected  to  be  a  linear  com¬ 
bination  of  derivatives  of  temperature  or  some  other 
pertinent  variable  of  interest,  i.e. 


or 


Calculations  of  the  flow  over  the  Viking  Aeroshell 
vehicle  using  this  scheme  are  presented  in  Ref.  13. 

Good  results  were  obtained  in  the  cases  reported  and 
sufficient  grid  adaption  was  used  to  adequately  resolve 
the  viscous  layer. 

The  adaptive  grid  calculations  cited  free*  Refs.  10 
and  13  are  for  adaption  along  either  a  constant  £  line 
or  a  constant  n  line.  These  one-dimensional  calcula¬ 
tions  can  be  achieved  using  either  Eq.  (23)  or  Eq.  (24) 
if  x  is  replaced  by  S  where  S  is  arc  lenqth  along  either 
constant  £  or  constant  0  lines.  In  addition,  S  can 
only  depend  upon  one  computational  coordinate.  For 
instance,  if  S  represents  arc  length  along  a  constant 
H  line,  then  S  is  assumed  to  only  be  a  function  of  £. 

In  this  case,  the  Dwyer  method  and  Gnoffo' s  approach  are 
identical.  However,  if  arc  length  along  £  or  D  surfaces 
in  physical  space  is  considered  to  depend  upon  both 
coordinates,  these  methods  are  not  the  same. 


In  particular,  suppose  the  analog  of  Eq.  (23)  along 
an  II  •  constant  surface  is  written 
S 
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If  a  derivative  is  taken  with  respect  to  S,  this  expres¬ 
sion  becomes 


(30) 


The  results  of  applying  Eq.  (23)  with  the  weighting 
function  of  Eq.  (25)  with  b  =  0  from  Ref.  10  is  shown 
in  Figs.  4  and  5.  In  this  problem,  the  equation  for 
unsteady  heat  conduction  was  solved  on  the  domain  shown 
with  the  temperature  set  equal  to  zero  everywhere .  ht 
t  -  0,  the  temperature  was  raised  impulsively  to  a  con¬ 
stant  value  on  the  lower  boundary  and  held  at  a  fixed 
value.  Figure  4  shows  the  isotherms  and  the  grid  at  an 
early  time  showing  the  high  temperature  gradients  near 
the  lower  boundary.  As  time  increases,  the  heat  flow 
into  the  domain  can  be  observed  in  both  the  isotherms 
and  the  grid  geometry  of  Fig.  5. 


where  the  right-hand  side  only  depends  upon  the  arc 
measure  along  the  constant  £  surfaces.  Differentiating 
again  yields 


The  identity  relating  £  and  the  usual  coordinate 
metrics  is 


ADAPTIVE  GRID  SCHEMES 


The  grid  adaption  is  provided  by  minimizing  a  weighted 
average  of  the  volume  variation  over  the  field  and  an 
appropriate  integral  is 


Brackbill  and  saltzman  (6.7)  and  Saltzman  and 
Brackbill  (23)  have  extended  Winslow's  method  (29)  for 
generating  a  computational  mesh  to  include  grid  adaption. 
In  this  approach,  a  variational  technique  is  used  to 
minimize  a  linear  combination  of  a  measure  of  grid 
smoothness,  orthogonality,  and  volume  (area)  variation. 
The  smoothness  is  swasured  by  integrating  the  change 
of  the  computational  coordinates  over  the  physical 
domain  and  for  two-dimensions  may  be  written 

I  -  /  [<V£)2  +  (Vn)2]dV  (7) 

*  D 

where  dv  represents  differential  volume  in  physical 
space  and  D  denotes  the  physical  domain.  The  smoothest 
mapping  between  the  physical  and  computational  domains 
is  obtained  by  minimizing  Ig  alone.  The  result  is 
that  Laplace's  equation  for  the  computational  coordi¬ 
nates,  £  and  g,  must  be  solved  to  determine  the  trans¬ 
formation. 

The  transformation  which  minimizes  IB  is  obtained 
by  solving  the  Euler-Lagrange  equations.  Equation  (7) 
is  first  written  as  an  integral  in  computational  space 
using  the  identities 

=  VJ  -  S  =  -VJ '  nx  *  "VJ  '  ny  "  VJ  (8) 


J  =  x.y  -  x  y, 

£  n  nr£ 


Equation  (7)  then  becomes 
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(10) 


where  indicates  integration  over  the  computational 
domain.  The  Euler-Lagrange  equations  for  this  varia¬ 
tional  problem  are 


If  the  indicated  differentiations  are  carried  out  in 
Eqs.  (11)  and  (12),  the  familiar  form  of  the  LaPlace 
grid  generator  in  the  computational  domain  is  obtained. 

In  addition  to  a  smoothness  requirement,  control 
of  mesh  skewness  is  obtained  if  a  measure  of  grid 
orthogonality  is  included.  This  orthogonality  Aeasure 
is  provided  by  the  integral 
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where  w  is  the  weighting  function  which  produces  the  grid 
adaption.  Clearly,  an  adaptive  system  is  obtained  if 
the  product  of  a  positive  weight  function  and  the  cell 
area  is  held  at  a  fixed  value  over  the  physical  domain. 


In  order  to  incorporate  smoothness,  orthogonality 
and  adaption  in  the  grid  generator,  the  linear  combina¬ 
tion  of  the  integrals  given  in  Eqs.  (10),  (14),  and  (16) 
is  written  as 
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If  the  Euler-Lagrange  equations  for  minimizing  Eq.  (17) 
are  derived,  they  are  of  the  form 
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where  the  aif  bi(  and  are  functions  of  the  weighting 
coefficients,!  and  X  .and  the  metric  coefficients  of 

'  V  O' 

the  transformation.  The  values  of  ai(  b^,  and  are 
given  in  Ref.  7  and  are  not  repeated  here.  The  mesh 
which  is  generated  when  Eqs.  (18)  are  solved  can  b>- 
varied  by  adjusting  the  weighting  coefficients.  For 
example,  a  large  value  of  Xy  will  provide  more  adaption 
in  the  grid  with  less  emphasis  on  orthogonality  and 
smoothness . 


Figures  2  and  3  show  the  pressure  contours  anti  the 
adaptive  grid  for  the  inviscid  supersonic  flow  over  a 
forward  facing  step.  The  similarity  between  the  pres¬ 
sure  contours  and  the  mesh  is  quite  good.  In  this  case, 
the  weight  function  for  grid  adaption  was  taken  to  be 
the  square  of  the  magnitude  of  the  pressure  gradient 
divided  by  the  pressure.  If  one  is  interested  in 
tracking  a  shock,  this  is  a  reasonable  choice.  Strictly 
speaking,  the  shock  wave  in  an  inviscid  flow  is  a  dis¬ 
continuity  and  even  with  grid  adaption,  the  shock 
mathematically  still  should  occur  between  two  grid 
points.  Of  course,  the  shock  is  given  a  pseudo-viscous 
profile  through  the  introduction  of  artificial  viscosity 
either  by  smoothing  the  data  or  by  the  form  of  the  algo¬ 
rithm.  In  this  case  a  flux-corrected  transport  scheme 
was  used.  This  provides  a  smooth  shock  profile  in  the 
solution. 
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The  extension  of  the  variational  approach  to  three- 
dimenaional  problems  is  straightforward.  Additional 
terms  must  be  added  to  include  the  third  dimension  but 
conceptionally  the  approach  is  the  same.  Saltzman  and 
Brackbill  (23)  have  presented  several  examples  of  three- 
dimensional  adaptive  grid  calculations  and  include  mesh 
generation  for  a  wing-fuselage  junction. 


another  good  reason  for  using  an  adaptive  grid.  Of 
course/  the  ideas  of  resolving  regions  of  rapid  change 
in  dependent  variable  and  reducing  error  are  not  mutually 
exclusive.  The  largest  numerical  errors  are  usually 
found  in  regions  where  the  solution  is  changing  most 
rapidly. 

The  idea  of  an  adaptive  grid  implies  th.  _  the 
solution  of  a  partial  differential  equation  is  being 
computed  using  some  sort  of  iterative  or  inarching 
technique  For  hyperbolic  or  parabolic  problems,  the 
solution  is  computed  by  advancing  in  space  or  time  and 
adjusting  the  grid  as  the  solution  progresses.  For 
elliptic  problems,  a  relaxation  procedure  provides 
intermediate  results  which  are  used  to  adjust  the  grid 
point  positions.  A  simple  linear  equation  can  be  used 
to  illustrate  some  of  the  considerations  that  must  be 
made  in  employing  a  solution  adaptive  grid. 

When  a  partial  differential  equation  is  transformed 
from  physical  to  computational  space,  the  metrics  of  the 
transformation  and  the  grid  speeds  appear  in  the  equa¬ 
tion.  These  additional  terms  must  be  evaluated  in  order 
to  solve  the  differential  equation  written  in  computa¬ 
tional  coordinates.  As  an  example,  the  firat-order  wave 
equation  may  be  written 

3u  3  u 

5t  +  c  sT  =  0  U) 

where  u(t,x)  is  the  unknown  dependent  variable  and  c  is 
a  constant  wave  speed.  For  this  one-dimensional  example, 
let  the  transformation  relating  the  physical  and  compu¬ 
tational  domains  be  written 


f,  -  f.(x.y) 

where  l  and  £,  are  coordinates  in  the  computational 
domain.  Transforming  the  wave  equation  into  computa¬ 
tional  coordinates  vields 

!?♦  <*t  +  eVaT-°  (3) 

The  terms  and  £  may  be  replaced  by  the  expressions 

*>t  =  _x/xf. 
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The  original  expression  may  then  be  written 


The  metric  represents  the  ratio  of  arc  lengths  in 
the  physical  and  computational  planes  and  the  grid  speed 
provides  the  dynamic  coupling  of  the  moving  grid  with 
the  evolving  solution  of  the  differential  equation. 

Any  method  for  constructing  an  adaptive  grid  must  pro¬ 
vide  e  means  of  estimating  these  terms  since  they 
explicitly  appear  in  the  transformed  differential  equa¬ 
tion.  Exceptions  to  this  can  be  cited.  In  problems 
where  a  time-asymptotic  solution  is  computed,  the  grid 
can  be  fixed  and  the  grid  speed  term  set  equal  to  zero. 
After  a  predetermined  number  of  iterations,  the  compu¬ 
tational  mesh  is  adjusted  and  the  iteration  process  is 
resumed.  This  regridding  procedure  is  equivalent  to 
solving  a  sequence  of  initial  boundary  value  problems. 
Other  exceptions  include  cases  where  the  solution 
changes  by  an  extremely  small  amount  over  a  single 
iteration  and  the  grid  speed  is  set  equal  to  zero  in 
the  governing  differential  equation.  However,  both  the 
grid  speed  and  metric  term  should  be  included  in  the 
general  case. 

Adaptive  grid  methods  can  be  divided  into  two 
categories.  In  the  first  category,  some  law  relating 
the  grid  points  in  the  physical  and  computational  domain 
is  used  to  establish  new  physical  grid  point  locations 
tx’s)  at  the  end  of  each  time  step.  The  grid  speed  term 
is  estimated  for  the  next  integration  step  by  using  a 
simple  backward  difference.  Hie  second  class  of  schemes 
relies  on  directly  establishing  the  grid  speed  by  some 
rule.  The  grid  speed  is  integrated  along  with  the  dif¬ 
ferential  equation  and  the  new  grid  point  positions  are 
established.  From  this  information,  the  metrics  are 
computed  by  evaluating  the  ratio  of  jure  lengths  in  the 
physical  and  computational  domains. 

Hiere  are  advantages  and  disadvantages  to  both 
approaches  for  generating  an  adaptive  grid.  Methods 
which  directly  generate  the  new  coordinates  through  a 
defined  mapping  are  conceptually  easy  to  apply.  Since 
the  grid  point  locations  are  established  through  the 
use  of  a  steady  grid  equation,  the  grid  speeds  are  most 
easily  determined  by  using  a  backward  difference.  This 
difference  is  usually  first  order  in  time  and  more 
accurate  dynamic  coupling  of  the  grid  motion  and  the 
partial  differential  may  be  desirable.  Some  grid  point 
location  schemes  initially  developed  for  one-dimensional 
applications  are  difficult  to  extend  to  two-  or  three- 
dimensional  problems.  Techniques  which  directly  deter¬ 
mine  grid  speed  from  some  grid  speed  law  are  easily  used 
in  multidimensional  applications  because  grid  point 
location  is  determined  by  a  simple  integration.  The 
major  disadvantage  is  that  physical  laws  which  relate 
grid  speed  and  grid  adaption  may  be  difficult  to  formu¬ 
late  and  the  success  of  the  method  depends  upon  the 
ingenuity  used  in  constructing  this  law.  For  these 
schemes,  point  control  is  also  a  problem. 


(XT  '  c)  3u 


When  a  solution  of  Eq.  (f>)  is  computed  using  a  fixed 


grid,  the  metric  coefficient. 


is  determined  from  the 


grid  geometry  which  is  established  at  the  beginning  of 
the  calculations.  This  term  does  not  change  bo  long  as 
a  fixed  mesh  is  used.  The  grid  speed  in  physical 
space,  x^  is  zero  in  this  case.  When  an  adaptive  grid 
is  used,  the  grid  speed  is  nonzero  because  the  grid 
changes  as  the  calculation  proceeds  and  the  metric 
coefficient  also  changes  each  time  the  grid  is  altered. 


A  number  of  adaptive  grid  generation  schemes  in 
both  classes  have  been  developed.  These  methods  have 
been  used  with  success  on  a  variety  of  problems.  In 
the  following  sections,  a  number  of  the  most  successful 
schemes  from  both  categories  are  reviewed.  In  this 
discussion,  the  problem  of  constructing  an  adaptive 
grid  it  viewed  as  one  of  allocation.  How  should  a  fixed 
number  of  grid  points  be  distributed  to  improve  the 
quality  of  a  numerical  solution?  This  distribution 
of  mesh  points  is  influenced  by  both  motion  of  the 
boundaries  and  solution  changes  on  the  interior  of  the 
physical  domain.  Hie  main  focus  of  interest  in  this 
paper  is  on  grid  point  motion  caused  by  solution  changes 
on  the  interior. 
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ABSTRACT 

A  number  of  techniques  of  constructing  adaptive 
mesh  generators  for  use  in  solving  partial  differential 
equations  are  reviewed  in  this  paper.  Techniques 
reviewed  include  methods  based  on  steady  grid  genera¬ 
tion  schemes  and  those  which  are  explicitly  designed 
to  determine  grid  speeds  in  a  time-dependent  or  space¬ 
marching  problem.  Results  for  candidate  methods  are 
included  and  suggestions  for  areas  of  future  research 
are  suggested. 


INTRODUCTION 

The  numerical  solution  of  the  partial  differential 
equations  governing  both  internal  and  external  flows  has 
reached  a  high  state  of  development  during  the  past 
fifteen  years.  Numerical  methods  for  solving  the  dif¬ 
ferent  types  of  equations  have  been  available  for  much 
longer.  However,  the  ability  to  treat  complex  geometries 
common  to  most  physical  problems  has  only  recently  been 
acquired.  In  fact,  the  pacing  item  in  advancing  numer¬ 
ical  procedures  for  solving  fluid  flow  and  heat  transfer 
problems  has  been  the  development  of  general  techniques 
for  numerically  constructing  mesh  systems  which  are 
boundary  conforming. 

The  construction  of  a  suitable  grid  is  the  first 
task  that  must  be  completed  when  the  numerical  solution 
of  a  system  of  partial  differential  equations  is  desired, 
dice  the  grid  is  generated,  the  system  of  equations  iB 
discretized  and  the  resulting  system  of  algebraic  equa¬ 
tions  is  solved.  This  solution  yields  the  values  of 
the  dependent  variables  at  each  of  the  mesh  points.  The 
solution  of  the  governing  equations  is  completed  in  a 
computational  domain  which  is  selected  to  be  rectangular 
shaped  for  simplicity.  The  physical  and  computational 
domains  are  related  through  a  mapping  as  schematically 
illustrated  in  Fig.  1  for  a  two-dimensional  case.  Th# 
problem  of  numerical  grid  generation  is  concerned  with 
techniques  for  establishing  this  relationship  between 
physical  and  computational  space.  Thompson  (24)  and 
Thompson  et  al.  (27)  have  presented  a  comprehensive 
review  of  the  state  of  the  art  in  numerical  grid  gener¬ 
al  ion . 

In  solving  partial  differential  equations  using 
numerical  methods,  the  selection  of  the  locations  of 
the  mesh  points  is  important  in  establishing  the  quality 
of  the  solution.  These  grid  point  positions  are  gener¬ 
ally  determined  initially  and  remain  fixed  throughout 
the  calculation.  In  order  to  determine  the  beat  grid 
point  locations,  an  a  priori  knowledge  of  the  solution 
of  the  physical  problem  is  desirable.  Unfortunately, 
this  knowledge  is  unavailable  and  only  the  general 
features  of  the  solution  may  be  initially  understood. 


For  example,  flow  over  a  body  must  be  computed  with  a 
grid  employing  a  sufficient  number  of  points  in  the 
viscous  regions  to  resolve  the  salient  features  of  the 
flow.  In  this  case,  a  grid  may  be  constructed  using  a 
compression  mapping  in  order  to  provide  a  large  number 
of  points  in  the  viscous  layer  near  the  body.  High 
mesh  densities  are  desirable  in  regions  where  large 
gradients  exist.  Since  the  exact  location  and  size  of 
these  regions  is  initially  unknown,  the  construction  of 
a  suitable  grid  in  the  general  case  is  difficult  and 
some  means  of  incorporating  information  from  the  solu¬ 
tion  in  locating  the  grid  points  is  needed. 

*nie  concept  of  a  solution  adaptive  grid  is  appealing 
for  a  number  of  reasons.  In  many  problems,  multiple 
length  scales  appear  and  a  grid  which  resolves  a  physical 
process  scaled  to  one  significant  length  can't  resolve 
events  which  occur  on  a  scale  less  than  the  size  of  the 
smallest  cell  or  mesh  increment.  A  typical  example  is 
again  provided  by  the  flow  of  a  viscous  fluid  over  a  body. 
The  inner  flow  near  the  body  in  the  boundary  layer  will 
be  resolved  in  sufficient  detail  to  be  of  use  only  if 
very  small  grid  spacing  is  used.  If  mesh  spacing  is 
used  with  the  minimum  size  determined  by  that  required 
to  resolve  the  outer  inviscid  flow,  the  detail  of  the 
boundary  layer  is  completely  lost.  Heat  transfer  and 
skin  friction  data  obtained  from  such  a  calculation  are 
cceipletely  meaningless.  With  the  use  of  an  adaptive 
grid,  the  physical  behavior  of  the  fluid  in  both  re- 
9lonB  be  adequately  established  using  the  same  set 
of  governing  equations.  The  different  length  scales  in 
the  problem  are  accoossodated  by  a  variable  mesh  size. 

In  a  sente,  this  approach  is  analogous  to  classical 
method*  which  require  a  solution  of  the  inner  and  outer 
flow  with  appropriate  matching  conditions.  Two  sets  of 
governing  equations  must  be  solved  while  numerically,  a 
single  set  of  governing  equations  is  solved,  but  the 
grid  position  problem  must  also  be  treated. 

When  a  parti.*’  differential  equation  is  discretized, 
errore  are  present  in  the  computed  solution.  If  the 
mesh  points  are  adjusted  during  the  calculation  to 
reduce  some  measure  of  the  local  or  global  error,  the 
quality  of  th#  solution  will  be  improved.  This  is 
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Fig.  17  Shock  aligned  grid  for  a  two-dimensional 
plane  shock  problem 
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Fig.  18  Comparison  of  confuted  pressure  solutions 
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Pig.  20  Grid  structure  and  corresponding  vorticity 
distribution  for  separated  flow,  R^  ■  100 
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Fig.  19  Grid  structure  and  corresponding  isotherm 
distribution  in  separated  flow,  R  -  100 
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ABSTRACT 

Methods  for  constructing  adaptive  grids  in 
more  than  one  dimension  have  been  developed. 

These  methods  are  usually  based  upon  the  idee  of 
equidistribution  of  a  weight  function  over  a  grid. 
Unfortunately,  for  large  grid  adaption  rates,  se¬ 
vere  skewing  occurs  in  the  mesh.  Two  techniques 
for  generating  an  orthogonal  adaptive  grid  are  de¬ 
veloped  and  results  of  applying  both  schemes  to 
some  simple  functional  examples  are  presented  for 
the  two-dimensional  case.  Extension  to  three  di¬ 
mensions  is  discussed  and  advantages  and  disadvan¬ 
tages  of  the  methods  are  identified. 


INTRODUCTION 

Grid  generation  has  always  been  a  problem  of 
major  concern  in  the  numerical  solution  of  partial 
differential  equations.  During  the  past  ten 
years,  satisfactory  methods  for  generating  body- 
t itted  mesh  systems  have  evolved  and  have  been 
used  with  great  success  on  a  variety  of  problems. 
More  recently,  a  great  amount  of  interest  has  cen¬ 
tered  on  the  development  of  dynamically  adaptive 
mesh  systems  which  evolve  with  the  solution  of  the 
PDE.  Adaptive  grid  schemes  are  attractive  and  are 
desirable  for  a  number  of  reasons.  These  reasons 
have  been  discussed  in  detail  by  a  number  of  au¬ 
thors  .  **  2 

Adaptive  schemes  in  one  dimension  have  been 
developed  and  applied  by  many  including  Gnoffo,3 
Dwyer  et  a  l.,*  and  Rai  and  Anderson.*  Basically, 
these  one -dimensional  schemes  all  rely  upon  equi¬ 
distribution  of  a  weight  function  over  a  mesh, 
i .  e.  , 


(1) 


where  w  is  some  positive  weighting  function,  is 

the  metric  of  the  transformation  from  physical  to 
computat ional  space,  and  C  is  a  constant.  Using 
this  law,  points  can  be  distributed  to  satisfy  any 
requirement  built  into  the  weight  function.  This 
expression  [Eq.  (1)1  was  solved  for  either  x  or  ( 
by  Dwyer  and  Gnoffo  using  a  direct  integration. 

Rai  and  Anderson  used  an  iterative  residual  ap¬ 
proach  to  obtain  the  same  result. 
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The  extension  of  the  equidistribution  idea  to 
more  than  one  dimension  is  desirable  since  most 
gains  from  dynamically  adaptive  grids  will  un¬ 
doubtedly  be  in  multidimensional  applications.  A 
logical  place  to  start  the  extension  to  two  dimen¬ 
sions  is  to  construct  a  two-coordinate,  indepen¬ 
dent  scheme  using  the  direct  integral  of  the  equi¬ 
distribution  law  given  in  Eq.  (1).  Anderson1  has 
reported  such  an  extension  and  has  applied  this  to 
simple  functions  to  study  the  adaption  process. 

The  main  difficulty  with  this  approach  is  that 
high  grid  skewness  occurs  even  for  moderate  grid 
adaption.  In  other  multidimensional  studies,  Rai 
and  Anderson  applied  their  scheme  to  a  number  of 
examples.  The  grid  adaption  employed  in  these  ex¬ 
amples  was  not  sufficiently  large  to  induce  skew¬ 
ness  problems. 

The  problem  of  controlling  grid  distortion  in 
constructing  adaptive  grids  for  multidimensional 
applications  must  be  addressed.  Methods  formulat¬ 
ed  without  a  direct  means  of  grid  skewness  control 
are  not  viable  in  applications  where  dense  point 
clustering  is  desired.  The  grid  distortion  prob¬ 
lem  is  avoided  if  grid  orthogonality  is  enforced 
when  a  mesh  is  generated.  In  this  paper,  two 
schemes  for  constructing  an  adaptive,  orthogonal 
mesh  are  presented.  While  these  methods  are  still 
in  the  exploratory /development  stage,  the  prelimi¬ 
nary  results  are  promising. 


PROBLEM  REVIEW 

In  order  to  understand  the  difficulty  of  con¬ 
structing  a  multidimensional  adaptive  grid,  it  is 
necessary  to  review  the  equidistribution  concept 
and  show  some  typical  results.  The  most  easily 
understood  multidimensional  scheme  employs  inde¬ 
pendent  grid  point  adaption  along  the  constant 
computational  coordinate  surfaces  in  physical 
space. 

Let  (t,n)  represent  the  computational  coordi¬ 
nates  and  (x,y)  be  coordinates  in  the  physical  do¬ 
main.  If  S  represents  arc  length  along  a  constant 
0  surface  in  physical  space,  a  simple  equidistri¬ 
bution  law  controlling  point  motion  along  this 
surface  may  be  written. 

S{wl  “  cl(,°  (2) 

where  Wj  is  a  positive  weight  function  which  de¬ 
pends  upon  the  solution  of  the  PDE  system  under 
consideration  and  c ^  is  a  function  of  the  computa¬ 
tional  coordinate,  n .  This  equation  can  be  solved 
to  obtain  the  integral  for  arc  length 


s/s 
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d£ 


ferential  equation  was  solved  with  a  numerical 
method  using  these  grids.  Some  means  of  control¬ 
ling  the  grid  skewness  must  be  incorporated  in  the 
mesh  generator. 


Method  0T1 


If  N  represents  art  length  along  the  £  equal  con¬ 
stant  surface,  tbe  companion  equation  with  weight 
function  is 


The  physical  coordinates  (x,y)  can  be  recov¬ 
ered  from  Eqs .  (3)  and  (4).  These  coordinates  can 
also  be  computed  directly  from  differential  equa¬ 
tions.  Since 


The  grid  skewness  at  each  point  in  a  mesh  is 
easily  evaluated  by  computing  the  angle  of  inter¬ 
section  between  constant  £  and  <)  surfaces.  Con¬ 
sider  the  Intersection  of  £  and  n  equal  constant 
lines  in  physical  space  (see  Fig.  3).  If  i^  is 

the  unit  vector  along  the  n  equal  constant  curve 
and  i^  is  the  unit  vector  along  the  £  equal  con¬ 
stant  curve,  the  cross  product  of  these  unit  vec¬ 
tors  may  be  written 


h "  V 


sin8 


where 
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(6) 


governing  differential  equations  for  (x,y)  can  be 
obtained  from  the  equid istribut ion  laws  along  £ 
and  n  equal  constant  surfaces.  For  example,  dif¬ 
ferentiating  Eq.  (2)  and  employing  Eq.  (5)  yields 
the  expression 
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(7) 


The  companion  expression  along  a  c 
face  may  be  written  (see  Ref.  2). 
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Typical  results  of  applying  this  independent 
equidistr ibution  concept  along  constant  $  and  r\ 
surfaces  are  shown  in  Figs.  1  and  2.  The  weight 
functions  were  of  the  form 


w,  -  1  +  A|g| 


W2  •  1  +  Bif,l 


(9) 


where  A  and  B  are  constants  which  determine  the 
magnitude  of  the  adaption  desired.  The  shock* like 
discontinuity  in  Fig.  1  is  typical  of  many  func¬ 
tions  where  rapid  changes  occur  along  one  primary 
coordinate.  The  choice  of  a  sinusoidal  shaped 
surface  in  Fig.  2  provides  functional  changes  in 
both  directions  and  leads  to  adaption  along  both 
families  of  coordinate  surfaces.  At  lower  values 
of  the  adaption  constants,  A  and  B,  distortion  is 
relatively  low.  However,  it  is  clear  that  signif¬ 
icant  distortion  occurs  in  both  cases  shown.  One 
would  not  expect  to  obtain  good  results  if  a  dif- 
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Performing  the  indicated  operations,  Eq.  (10)  may 
be  written 

J  =  S^sinB  (12) 

where  6  is  the  intersection  angle  between  £  and  n 
equal  constant  curves  snd  J  is  the  Jacobian  of  the 
trans format ion 

J  *  x£yn  '  V{  (13) 

The  intersection  angle  is  easily  monitored  by  com¬ 
puting  sin8  through  Eq.  (12).  In  fact,  the  angle 
8  is  directly  influenced  by  the  choice  of  weight 
functions.  If  the  equid istribut ion  laws  for 

and  are  substituted  into  Eq.  (12),  the  result 

is 

c1c2/w,w2  =  J/sinB  (14) 


Assume  that  grid  adaption  along  one  coordinate 
(S)  is  all  that  is  necessary.  Since  the  arc  ele¬ 
ments  along  the  S  direction  are  calculated  inde¬ 
pendently,  the  arc  is  given  by  Eq.  (3).  Suppose 
that  grid  distortion  is  controlled  by  selecting 
the  intersection  angle,  8  in  Eq.  (14).  However  if 
8  is  specified,  the  constant  c^  and  the  weight 

function  Wj  cannot  be  independently  chosen.  Since 
Cj/Wj  and  8  ara  given, 


c2/w2 


J(w1/c1) 

sinS 
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and 


ll.  *  i  I  =  sin6 
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(19) 
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(16) 


or 


|(lx{  +  jy£  +  kz£)  x  (ix^  +  Jyn  +  kz^l 


The  constant  is  of  the  form 


It  is  interesting  to  note  that  the  length  scale 
provided  through  does  not  appear  in  Eq.  (16). 

This  shows  that  the  value  of  N  computed  at  a  given 
point  is  not  scaled  as  in  £q.  (4).  The  absence  of 
a  normalizing  length  scale  in  Eq.  (16)  indicates 
that  the  system  of  PDEs  governing  this  scheme  is 
hyperbolic.  The  solution  of  such  a  problem  must 
be  computed  in  the  computational  domain  by  speci¬ 
fying  initial  data  at  i)  =  0  (N  =  0),  and  marching 

the  solution  outward  to  n  The  outer  boundary 

max 

in  physical  space  corresponding  to  n  must  float 

max 

and  is  determined  as  part  of  the  solution. 

Figure  4  shows  a  solution  for  the  same  shock¬ 
like  function  employed  in  Fig.  1  in  the  physical 
domain.  In  this  case,  the  angle  8,  has  been  se¬ 
lected  to  be  90  degrees  so  an  orthogonal  mesh  is 
created.  The  grid  is  21  *  21  and  the  weight  func¬ 
tion  Wj  is  given  in  Eq.  (9).  Again  the  outer 


=  S,N  sin8  (20) 

The  angle  between  the  normal  to  the  plane  formed 

by  the  unit  vectors  i.  and  i  and  the  unit  vector 
.  4  9 

can  be  monitored  by  forming  the  box  product 

(i£  x  i^)  .  =  sin*  (21) 

This  expression  reduces  to 

StN  M  =  J/sin8  (22) 

4  9  4 

If  it  is  assumed  that  the  adaption  is  in  the  (  di¬ 
rection,  the  quantity  S£  is  prescribed  by  the 

equidistribution  law.  Equations  (20)  and  (22) 
provide  the  additional  expressions  for  the  quanti¬ 
ties  N  and  M, .  This  formulation  is  similar  to 
9  4 

the  two-dimensional  case.  While  the  linearization 
and  classification  of  the  system  for  three  dimen¬ 
sions  has  not  been  done,  it  is  reasonable  to  ex¬ 
pect  that  the  grid  would  be  computed  on  an  open 
domain.  The  governing  equations  are  probably  hy¬ 
perbolic  due  to  the  similarity  to  the  two-dimen¬ 
sional  case.  An  orthogonal  grid  can  be  generated 
with  adaption  in  one  direction  if  *  and  6  are  both 
taken  to  be  90  degrees. 


boundary  is  free  and  the  solution  determines  the 
final  shape.  Both  the  orthogonality  and  the 
floating  outer  boundary  are  apparent  in  these 
results. 

Figure  5  shows  the  orthogonal  grid  generated 
using  the  sinusoidal  function  of  Fig.  2.  The  re¬ 
sulting  grid  shows  no  evidence  of  skewing  although 
the  distorted  outer  boundary  is  again  apparent. 

In  almost  all  calculations  involving  dynamically 
adapting  mesh  systems,  the  grid  point  speeds,  lo¬ 
cations,  or  the  forcing  functions  are  smoothed. 

In  computing  adaptive  grids  with  0T1,  it  was  Ob¬ 
served  that  adding  smoothing  relaxed  the  orthogo¬ 
nality  condition.  Thus,  no  smoothing  was  employed 
in  computing  the  results  for  method  0T1. 

The  results  for  the  orthogonal  calculations 
shown  are  very  good.  In  applications  where  adap¬ 
tion  is  one  dimension  is  desirable  and  a  free  out¬ 
er  boundary  is  not  a  problem,  this  is  a  viable  ap¬ 
proach  for  generating  an  adaptive  grid.  Can  this 
scheme  be  extended  to  three  dimensions? 

In  three  dimensions,  the  applicable  equidis- 
trihuion  laws  would  be  of  the  form 


h  =  ci/wi 

(18a) 

N!)  =  C2/W2 

(18b) 

=  C3/W3 

(18c) 

where  S,  N,  and  M  are  arc  lengths  along  the  compu¬ 
tational  coordinate  surfaces.  The  angle  between 
the  4  and  e  directions  can  again  be  controlled  by 
noting  that 


Hethod  0T2 

The  method  presented  above  is  most  attractive 
for  generating  adaptive  grids  for  those  problems 
where  grid  adaption  is  necessary  in  only  one  di¬ 
rection.  However,  it  seems  more  appropriate  to 
employ  an  equidistribution  law  based  on  cell  area 
or  volume  rather  than  arc  length  when  problems  in 
more  than  one  dimension  are  considered.  In  a  re¬ 
cent  paper,  Anderson*  has  introduced  such  an  idea 
where  the  equidistribution  law  is 

Jw  *  Ap(t)/Ac  (23) 

In  this  expression,  w  is  a  positive  measure  of  the 
solution  and  is  the  weight  function,  J  is  the  Ja¬ 
cobian,  Ap(t)  is  the  physical  domain  integral 

Ap(t)  *  H  wdxdy  (24) 

and  Ac  is  the  area  of  the  computational  domain. 

This  equidistribution  law  is  incorporated  in  an 
area  continuity  equation 


where  (xT>  y^)  are  the  grid  speeds.  If  an  orthog¬ 
onal  grid  is  desired,  the  time  derivative  of  the 
orthogonal  condition 

a/»T  [x^  +  y£y^ )  =  0  (26) 


3 


SUMMARY  AND  CONCLUSIONS 


completes  the  set  for  the  unknowns  (x^ ,  y^).  ^ 

the  initial  grid  is  orthogonal,  the  final  grid 
will  be  orthogonal.  After  the  grid  speeds  have 
been  determined,  the  grid  point  locations  are  ob¬ 
tained  by  integration  with  respect  for  T. 

The  system  of  Eqs .  (25)  and  (26)  is  hyperbolic 
and  the  grid  speeds  are  determined  by  marching  the 
solution  outward  away  from  the  initial  data  sur¬ 
face  at  n  =0.  Boundary  conditions  can  be  en¬ 
forced  at  £  =  0  and  £  =  £  The  outer  boundary 

max 

is  free  to  float  as  determined  by  the  integration 
of  the  grid  speeds.  Notice  that  this  system  is 
weakly  elliptic  through  the  source  term  of  Eq. 

(25).  For  the  results  shown  in  Figs.  6  and  7,  the 
outer  boundary  was  held  at  a  fixed  position  and 
the  last  grid  line  computed  in  the  hyperbolic 

marching  scheme  was  at  t)  -  Art  .  The  source  term 
®  max 

was  computed  on  a  domain  with  fixed  boundaries. 

The  results  show  that  this  grid  is  also  adaptive 
and  orthogonal.  In  this  case  the  weight  function 
was  selected  to  be  of  the  form 

w  =  1  +  | Vr  u |A  (27) 

£  »n 

The  results  in  Figs.  6  and  7  for  the  shock  and 
sine  function  are  computed  with  grid  adaption 
based  upon  an  area  equidistr ibut ion .  While  some 
similarities  exist  between  the  results  obtained 
using  0T1  and  0T2 ,  the  grids  produced  do  show  some 
differences.  One  of  the  problems  noted  when  em¬ 
ploying  0T2  is  that  grid  adaption  is  always  accom¬ 
plished  at  the  expense  of  available  cell  area  at  a 
greater  value  of  T| .  Since  the  grid  equations  are 
hyperbolic,  the  area  equidistr ibut ion  law  always 
necessitates  the  borrowing  of  area  at  large  n. 
Consequently,  the  adjustment  of  the  mesh  is  slow 
since  the  source  term  of  Eq .  (25)  is  the  only 
means  for  providing  an  upstream  influence. 

For  two-dimensional  problems,  the  grid  pro¬ 
duced  by  either  method  0T1  or  0T2  are  satisfactory 
for  the  cases  considered.  The  extension  of  0T2  to 
three  dimensions  can  be  accomplished  in  a 
straightforward  manner.  The  orthogonality  condi¬ 
tion,  Eq .  (26),  must  be  altered  to  include  the 
term  z_z  .  In  addition,  another  condition  is  nec- 

£  n 

essary  to  provide  the  proper  cell  orientation. 

This  relationship  is  supplied  by 

3/  •  r  (J  -  S  N  M  )  =  0  (28) 

£  n  C 

This  last  expression  provides  another  PDE  for  the 
grid  speeds  in  three  dimensions. 

ft  is  interesting  to  note  that  one  can  derive 
the  relationship  between  the  area  weight  function, 
w,  and  the  one-dimensional  functions  w^  and  w^ . 

For  the  orthogonal  grids  considered,  the  two-di¬ 
mensional  result  is 

h\  =  j  =  V2/W1W2  =  Vt)/wAc  (29) 

Since  ci'*;  ls  selected  and  c^/w0  Is  determined  by 

orthogonality,  the  corresponding  w  for  the  area 
equidistribution  case  can  be  directly  computed 
from  Kq  .  (29  ) . 


Two  schemes  for  producing  adaptive,  orthogonal 
grids  have  been  presented.  The  first  is  based 
upon  one-dimensional  equidistribution  and  provides 
adaption  along  only  one  coordinate.  The  other  co¬ 
ordinate  location  is  determined  by  the  orthogonal¬ 
ity  constraint.  The  second  method  employs  the 
concept  of  equidistribution  over  an  area  or  volume 
to  generate  a  single  PDE  for  the  grid  speeds  or 
point  locations.  Additional  expressions  are  ob¬ 
tained  from  the  orthogonality  conditions. 

Both  schemes  produce  systems  of  hyperbolic 
partial  differential  equations.  This  is  expected 
since  even  in  the  general,  nonadaptive  case,  or¬ 
thogonal  grids  cannot  be  obtained  on  a  closed  do¬ 
main  when  Dirichlet  boundary  conditions  are  used 
in  solving  the  governing  PDEs. 

Method  0T1  can  be  implemented  by  solving  for 
arc  lengths  along  constant  £  and  n  surfaces  and 
then  computing  the  corresponding  values  of  x  and 
y.  However,  an  alternative  approach  is  to  compute 
x  and  y  directly  from  the  governing  PDEs. 

Method  0T2  was  formulated  using  grid  speeds. 
However,  a  steady  formulation  may  also  be  used. 
With  this  approach,  the  governing  linearized  PDEs 
must  be  solved  by  marching  outward  away  from  an 
initial  data  surface.  In  thiB  case,  the  x  and  y 
coordinates  are  obtained  instead  of  the  grid 
speeds.  The  steady  formulation  of  method  0T2  is 
exactly  the  adaptive  counterpart  of  the  grid  gen¬ 
eration  scheme  presented  by  Steger  and  Sorenson.’ 

Both  schemes  provide  reasonable  results  for 
the  simple  test  problems  illustrated  here.  Addi¬ 
tional  work  ls  needed  in  evaluating  the  applica¬ 
bility  of  these  ideas  to  actual  flow  problems. 
Studies  coupling  the  orthogonal  generators  with 
the  flow  equations  will  commence  in  the  near  fu¬ 
ture. 
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Figure  2.  Grid  adaption  for  sinusoidal  function 
distribution,  A»  2,  B  -  3 

u(x.y)  -  0  0  <  y  <  9.  +  ltain(~) 

u(x,y)  -  0. 5 Jy  -  9.  -  4sin(~)  J  . 

9.  +  4sin (y-p  )  <  y  <  11.  +  4ein(-yy) 
u(x,y)  -  1.0  11.  +  4ain  yy  <  y  <  20. 


Figure  4.  Orthogonal,  one-dimensional  adaptive  grid 
for  shock-like  function  (see  Fig.  1), 

A  -  3.0 
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Figure  5. 


Orthogonal,  one-dimensional  adaptive  grid 
for  sinusoidal  function  (see  Fig.  2), 

A  -  3.0 


Figure  6.  Orthogonal,  two-dimensional  adaptive  grid 
for  shock-like  function  (see  Fig.  1), 

A  -  1.0 


Figure  7.  Orthogonal,  two-dimensional  adaptive  grid 
for  sinuaoidal  function  (see  Fig.  2), 

A  -  1.0 
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ABSTRACT 

An  approach  to  solution  adaptive  grid  genera¬ 
tion  for  use  with  finite  difference  techniques, 
previously  demonstrated  on  model  problems  in  one 
space  dimension,  has  been  extended  to  multidimen¬ 
sional  problems.  The  method  is  based  on  the  popu¬ 
lar  elliptic  steady  grid  generators,  hut  is  "dy¬ 
namically”  adaptive  in  the  sense  that  a  grid  is 
maintained  at  all  times  satisfying  the  steady  grid 
law  driven  by  a  solution-dependent  source  term. 

Testing  has  been  carried  out  on  Burgers'  equa¬ 
tion  in  one  and  two  space  dimensions.  Results  ap¬ 
pear  encouraging  both  for  inviscid  wave  propaga¬ 
tion  cases  and  viscous  boundary  layer  cases, 
suggesting  that  application  to  practical  flow 
problems  is  now  possible.  In  the.  course  of  the 
work,  obstacles  relating  to  grid  correction, 
smoothing  of  the  solution,  and  elliptic  equation 
solvers  have  been  largely  overcome.  Concern  re¬ 
mains,  however,  about  grid  skewness,  boundary  lay¬ 
er  resolution  and  the  need  for  implicit  integra¬ 
tion  methods.  Also,  the  method  in  3-1)  is  expected 
to  be  very  demanding  of  computer  resources. 


NOMENCLATURE 

a  =  clustering  constant 

c,d  =  wave  speeds 

e  =  numerical  error 

f.g  -  flux  vectors  in  governing  equation 

F  =  function  describing  surface  (F=0) 

G  =  grid  speeds  =  (x  .y^) 

vJ  =  Jacobian  of  trans format  ion 

k  -  the  quantity  J(Pr*  +  Or  ) 

s  n 

P.Q  -  grid  clustering  (forcing)  functions 

f  ~  position  vector  -  (x,y) 

R  =  residual  of  steady  grid  equation 

s  -  distance  along  surface 

t  =  l i me 

u  =  model  dependent  flow  variable 

w  -  positive  measure  of  solution  gradient 

(x,y)  -  physical  coordinates 

a, &,2l  -  coefficients  arising  f rom  transformation 

of  Laplacian 

-  computational  coordinates 
*  =  damping  constant 

P  -  viscosity  coefficient 

t  •  « .ompuir.t  iona  I  time 

w  -  smoothing  constant 
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V  =  backward  difference  or  gradient  operator 

A  =  forward  difference  or  increment  (as  in 

2  AS,  Ax) 

=  Laplacian  operator  in  x,y  domain 
=  Laplacian  operator  in  S.n  domain 

Matrices 

(A]  =  coefficient  matrix  in  grid  speed  equa¬ 

tion 

(B)  =  coefficient  matrix  in  grid  speed  equa- 

t  ion 

[Pu]  =  Jacobian  matrix  of  derivatives  of  P  with 
respect  to  u 

[Qu]  =  Jacobian  matrix  of  derivatives  of  Q  with 
respect  to  u 

( S ]  =  Smoothing  operator  expressed  in  matrix 

form 

Subscripts 

i,j  =  row  and  column  indices 

k,l  =  summation  indices 

x,y,z  =  partial  differentiation 

S.n.T  =  partial  differentiation 

Superscripts 

~  =  smoothed  quantities 

=  simplified  forms 


INTRODUCTION 

Methods  for  generating  fixed  finite  difference 
grids  around  two-dimensional  airfoils  and  other 
geometries  have  evolved  to  the  point  where  such 
grids  are  routinely  employed.  Often  these  grids 
are  generated  by  an  elliptic  partial  differential 
equation  relating  the  physical  and  the  computa¬ 
tional  domains'  (see  Chapter  10).  The  need  for 
developing  a  class  of  solution-adaptive  methods 
may  arise  from:  (1)  boundary  motion  in  unsteady 
flow  problems;  (2)  moving  shock  problems;  (3) 
time-marching  to  a  steady  state,  where  regions  re¬ 
quiring  high  resolution  are  not  known  in  advance; 
(4)  space-marching  problems  (e.g.  parabolized  Nav- 
ier-Stokes) ,  where  a  grid  must  be  generated  in 
each  transverse  plane  moving  downstream. 

It  is  natural  to  suppose  that  the  steady 
(static)  grid  generators  might  be  extended  to  the 
dynamic  case.  This  can  be  accomplished  by  differ¬ 
entiating  the  elliptic  p.d.e.  with  respect  to  time 
to  yield  the  so-called  grid  speed  equation  de¬ 
scribing  point  motion.  Grid  speeds  are  then  inte¬ 
grated  in  time  simultaneously  with  the  solution  to 
the  governing  equation  at  the  new  locations.  Such 
an  approach  was  tested  by  Hindman,  Kutler,  and  An¬ 
derson  on  an  Euler-equation  solver  in  two  dimen¬ 
sions  for  the  case  of  arbitrary  boundaries  but  no 
interior  grid  clustering,  and  more  recently  by 
Hindman  and  Spencer’  on  Burgers'  equation  in  one 
dimension  with  a  source  term  for  clustering. 


The  present  work  is  a  refinement  of  Hindman's 
1-D  method,  followed  by  an  extension  to  multidi¬ 
mensional  problems.  It  should  be  stressed  that 
the  main  focus  of  this  study  has  been  on  method 
development  rather  than  application,  hence  only 
simple  geometries  have  been  considered.  As  will 
be  discussed  later,  application  to  practical  1-D 
and  2 -D  problems  now  appears  possible. 

Headers  who  are  interested  in  a  more  complete 
introduction  to  adaptive  grid  methods,  or  who  wish 
to  study  alternate  approaches,  should  consult  the 
recent  survey  papers  by  Anderson . *  ~7 


METHOD 

In  order  to  efficiently  summarize  the  method, 
a  table  of  equations  is  provided  (Table  1). 

Please,  refer  to  this  table  while  reviewing  the 
following  comments.  Also  refer  to  the  block  dia¬ 
gram  of  the  "system*’  formed  by  the  grid  and  the 
model  flow  variable  (Fig.  1). 


Doma i n 

The  generalized  transformation  mapping  the 
physical  space  into  the  computational  space  is  as 
given  by  Eq  set  (1).  The  purpose  of  the  trans¬ 
formation  is  to  allow  a  uniform  red  angular  grid 
to  be  used  for  computation,  while  the  physical 
grid  conforms  to  boundary  shapes  and  is  clustered 
where  high  resolution  is  needed.  As  a  result  of 
the  mapping,  derivatives  must  be  transformed  ac¬ 
cording  to  Eq.  set  (2). 


Steady  Grid  Law 

A  specific  form  of  the  above  mentioned  trans¬ 
formation  is  obtained  by  solving  a  Poisson  equa¬ 
tion  (Eq.  3a),  an  idea  developed  by  Winslow  and 
Thompson9'1  °f or  the  t  ime- invariant  case.  To  solve 
for  the  (x,y)  point  locations,  the  role  of  depen¬ 
dent  and  independent  variables  must  be  inter¬ 
changed,  yielding  in  2-1)  a  coupled  set  of  two  non¬ 
linear  elliptic  partial  differential  equations 
(Eq.  3b).  The  forcing  I  unctions  P(x,y,u)  and 
Q(x,y,u)  are  carefully  chosen  to  provide  adequate 
clustering  without  grid  crossing  or  overlap. 


Cluster ing  Function 

A  rationale  for  choosing  P(x,u)  in  1-D  is  pro¬ 
vided  by  the  integral  grid  law  cited  (Eq.  4), 
where  the  function  w(u)  is  a  positive  measure  of 
solution  gradient  Differentiation  twice  with  re¬ 
spect  to  £  reveals  that  this  form  is  equivalent  to 
the  P-function  selected  (Eq.  S).  It  is  clear  that 
x  will  always  be  a  mono ton ica 1 ly  increasing  func¬ 
tion  of  £,  hence  grid  crossing  cannot  occur  (ex¬ 
cept  due  to  numerical  effects).  The  clustering 
constant  controls  the  amount  of  adaption,  from 
none  (a  =  0)  up  to  i  large  amount  (a  -->  ••)  .  The 
possibility  of  using  a  different  clustering  func¬ 
tion  will  be  discussed  later  In  2-1),  P  and  Q  are 
obtained  by  a  logi<al  extension  of  the  1-D  form. 


Simplified  forms  of  the  grid  equations  (Eqs. 

6,  7b)  resulting  from  the  cancellation  of  (x^) 

factors  were  used  in  1-D.  As  a  result,  the  equa¬ 
tions  are  linearized  and  may  be  solved  directly. 
In  2-D,  P  and  Q  can  still  be  split  into  a  factor 
involving  only  u  and  one  involving  only  x  and  y. 


Grid  Speed  Equation 

This  equation  (Eq.  7a, b)  describing  the  rates 
of  point  motion  (x^,y^)  is  derived  by  differenti¬ 
ating  the  steady  grid  law  with  respect  to  time 
(t).  Grid  speeds  would  be  zero  (with  no  boundary 
motion)  except  for  the  fact  that  changes  in  the 
solution  cause  the  clustering  functions  to  vary. 
Thus  it  is  very  important  that  accurate  represen¬ 
tations  for  and  be  obtained.  In  the  present 

work,  this  is  done  by  taking  analytical  deriva¬ 
tives  (Table  3).  Usually  P  and  Q  depend  only  on 
x,y,  and  u  at  the  center  point  of  the  finite  dif¬ 
ference  molecule  and  its  immediate  neighbors, 
which  simplifies  the  calculation  of  and  . 

Even  so,  substantial  amounts  of  computation  and 
storage  are  required,  suggesting  that  it  may  be 
better  to  obtain  these  quantities  by  backward  dif¬ 
ferencing  in  time  as  part  of  an  iterative  process 
for  solving  the  grid  speed  equation. 


Governing  Equation 

The  p.d.e.  governing  the  physical  process  is 
cast  into  conservation  law  form  (Eq.  8a),  then 
transformed  to  computat ional  space,  including 
terms  due  to  grid  point  motion  (Eq.  8b).  For 
Burgers'  equation  (Eqs.  8c, d)  the  flux  vectors  f 
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and  g  take  the  scalar  form  f  =  g  =  u  /2.  Linear¬ 
ized  expressions  f  =  cu  and  g  =  du  were  used  for 
simplicity  in  2-D.  In  the  case  of  the  Navier- 
Stokes  or  Euler  equations,  u,  f,  and  g  form  a  set 
of  vector  quantities,  but  only  one  element  of  u 
(say  the  density)  need  be  selected  to  drive  the 
grid.  All  flux  derivative  terms  for  the  inviscid 
problem  were  evaluated  by  upwind  differencing, 
analogous  to  the  flux  splitting  class  of  methods 
applied  to  the  Euler  equations. 

A  special  problem  is  presented  by  the  viscous 
Burgers'  equation  in  that  upwind  differencing  on 
the  convective  terms  generates  excessive  dissipa¬ 
tion  and  hence  large  errors  in  the  steady-state 
solution.  This  may  be  overcome  by  using  a  cen¬ 
trally-differenced  scheme,  however  grid  size  is 
then  restriced  to  satisfy  a  mesh  Reynolds  number 
constraint.  The  third-order  upwind  correction 
proposed  by  Leonard11  was  found  to  be  successful 
in  one  dimension.  (Also  see  results  section.) 


ALGORITHM 

The  computer  program  written  for  this  study 
executes  the  steps  that  follow  with  the  explana¬ 
tion  tailored  for  the  2-1)  case.  The  physical  do¬ 
main  in  2-D  is  the  unit  square. 
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Smoothing 


1.  Establish  input  data:  (a)  initial  distribution 
of  u  as  a  function  of  (x,y)  or  (£,n);  (b)  cluster¬ 
ing  constant,  smoothing  constant,  and  grid  damping 
constant  (see  below  for  explanation  of  smoothing 
and  damping);  (c)  time  step  limits  and  wave  speeds 
if  applicable;  (d)  iteration  limits,  tolerances, 
and  over-relaxation  factors;  (e)  number  of  steps 
to  be  computed  and/or  steady-state  convergence 
criterion;  (f)  initialization  of  (x,y)  to  uniform 
grid,  and  and  to  zero. 

2.  Solve  steady  grid  generator  for  initial  grid. 

A  Gauss-Seidel  point  iterative  process  is  em¬ 
ployed,  updating  x,y,P,Q  and  related  quantities 
during  each  sweep  through  the  grid.  Over-relaxa¬ 
tion  is  generally  possible  and  speeds  convergence 
considerably.  Boundaries  of  the  2-1)  domain  are 
treated  by  applying  the  1 -I)  steady  grid  law.  The 
forcing  functions  P  and  Q  may  be  calculated  in  ad¬ 
vance  if  u  is  initially  specified  as  a  function  of 

otherwise  ulx.y)  must  be  corrected  by  in¬ 
terpolation  or  other  means,  and  P  and  Q  recalcu¬ 
lated  as  (x,y)  change  at  a  point. 

Begin  loop 

3.  Cali. u late  the  transform.it  ion  metrics 

(x_,x  ,y.,v  ),  the  coefficients  a,  0,  J,  and  the 
S  i  t  n 

Jacobian  J.  Also  calculate  the  coefficient  matri¬ 
ces  A  and  H  in  the  grid  speed  equation  (Table  2) 

and  the  residuals  KX  and  of  the  steady  grid 
genet  dot  . 

4.  Solve  Liu*  gr  id  speed  equation:  (a)  calculate 

the  six  nonzero  P  and  ()  quantities  at  each 
u  it  ' 

point  ,  h«‘g  *  1 1  loop  U) )  calculate  ti^  ,  which  depends 
upon  C  x t  . y ^  ) ,  the  t  lux  vectors,  and  the  viscous 

trims  if  piesent;  (cl  update  the  boundary  grid 
spends  by  a  line  relaxation  scheme  acting  on  the 
l-l>  grid  speed  equation;  (d)  update  the  interior 
grid  speeds  by  an  alternating  di lection  line  re  - 
I  Uxrtt  ;on  scheme,  whereby  sweeps  in  the  ^-direction 
update  t  lie  x-equat  ion  and  sweeps  in  the 
n-di  i  »'c  i  ton  update  the  y -equal  ion  ;  (e)  exit  loop 

i  1  .  f  i.’  e  r  geil 

Notes'  l  tide  r  - 1  e  1  axat  lot,  must  he  used  to 
achieve  converged* c ,  especially  f  j  large  cluster¬ 
ing  constants.  It  is  important  to  include  the  de¬ 
pend*  in  *  <*j  on  (xt  ,y  l  implicitly  in  the  relax¬ 
ation  Irnir  w  1 1  *  ■  r  e  V  e  i  pn.sibl**. 

3  'ibtiin  m*w  flow  solut  ion  and  grid  lot. at  ions: 
(a)  re  -  *  a  I  <  a  I  ate  for  converged  gri*i  speeds;  (b) 

establish  allowable  time  step  size  » limited  by  an 
tnvtscid  GP1,- type  condition  and  by  a  viscous  term 
condition);  (c)  perform  first-order  explicit  inte¬ 
gration  to  get  new  values  tor  x,y,  and  u;  (d)  up¬ 
date  the  functions  i  and  ■)  at  the  new  time  level. 

Exit  loop  if  requested  number  of  tune  steps  have 
been  completed,  or  if  solution  has  reached  a 
steady  state. 


There  are  at  least  three  reasons  why  smoothing 
of  the  solution,  u,  might  be  necessary:  (a)  to 
suppress  oscillations  ("wiggles")  associated  with 
certain  methods  of  integrating  the  governing  equa¬ 
tion;  (b)  to  avoid  grid  crossing  that  may  occur 
due  to  numerical  effects  when  the  forcing  func¬ 
tions  become  large;  (c)  to  allow  the  derivatives 
and  to  be  computed  analytical ly ,  even  near  a 

discontinuity  in  u.  For  the  present  algorithm 
only  the  last  two  reasons  apply,  since  integration 
schemes  are  selected  to  avoid  oscillations  in  u. 
Therefore,  smoothing  is  needed  only  for  the  pur¬ 
pose  of  grid  calculations,  an  important  distinc¬ 
tion  because  the  undesirable  effects  of  artificial 
smoothing- induced  diffusion  and  dispersion  are 
avoided  when  integrating  the  governing  equation. 


To  understand  the  nature  of  problem  (b),  con¬ 
sider  the  simplified  grid  equation,  xu  +  xsp(u)  = 

0.  It  is  easily  seen  that  discretization  by  cen¬ 
tral  differencing  on  x^  and  x^  will  cause  grid 

crossing  if  abs(P)  >  2,  even  though  the  exact 
mathematical  solution  does  not  exhibit  this  prop¬ 
erty.  However,  by  applying  a  smoothing  operator, 
one  can  always  prevent  the  forcing  function  from 
becoming  too  large  and  also  insure  that  sufficient 
smoothness  exist  to  compute  and  0^ . 


Table  4  presents  1-D  and  2-D  versions  of  the 
smoothing  operator.  Since  both  (P^J  and  [S]  are 

tridiagonal  in  1-D,  the  grid  speed  equation  (con¬ 
taining  )  becomes  pentadiagonal ,  requiring  2.5 


times  as  much  computation  to  solve  as  the  standard 
Thomas  algorithm  for  tridiagonal  systems.  In  2-D, 
the  equations  are  solved  iteratively  and  one  need 
only  smooth  u^  before  updating  the  grid  speeds  on 

each  step.  For  the  viscous  Burgers'  equation,  no 
smoothing  is  required  if  the  starting  solution  is 
smooth.  In  practice,  such  an  initial  condition 


may  be  obtained  by  solving  V  =  0,  provided 

that  one  is  not  interested  in  transient  behaviors 
associated  with  other  possible  initial  conditions. 


Grid  Correction 

The  grid  speed  equation  may  bo  written  in  the 
form  dK/dt  =  0,  whore  R  is  the  residual  of  the 
steady  grid  generator.  Such  a  form  is  neutrally 
stable,  that  is  errors  in  integration  will  neither 
be  amplified  nor  damped  out.  An  obvious  fix  is  to 
add  a  damping  term,  dR/dT  +  XR  =  0,  where  the 
damping  constant  X  may  be  chosen  as  large  as  1/At 
for  stability  with  an  explicit  integration  scheme. 
The  resulting  grid  correction  method  is  both  ef¬ 
fective  and  much  simpler  than  previous  methods 
(involving  a  new  solution  to  the  steady  grid  equa¬ 
tion  while  holding  u  fixed  in  either  physical  or 
computational  space). 


SPECIAL  CONCERNS 

Although  tin-  algorithm  is  quite  straightfor¬ 
ward.  a  tew  special  points  need  to  be  addressed: 


Boundary  Conditions 


Dirichlet  boundary  data  can  be  specified  by 
fixing  (x,y)  at  points  along  a  £=const  or  i)=const 
boundary.  Unfortunately,  this  does  not  allow  for 


grid  adaptation  on  the  boundary.  A  better  idea  is 
to  use  the  I -D  grid  method  to  treat  the  bounda¬ 
ries.  busy  application  of  the  1 -I)  method  is  pos¬ 
sible  lor  st  r  .light  boundaries  (as  considered 
here);  however,  for  curved  boundaries,  a  more  gen¬ 
eral  form  must  be  used,  +  s^P(u)  =  0,  where  s 

is  ill**  distance  along  the  surface.  In  the  grid 
spc“d  equation  the  general  boundary  condition  is 
G  •  VK  =  ~3F/3t,  where  G  is  the  vector  (x^.y^)  of 

grid  speeds,  and  P(x,y,t)  =  0  describes  a  moving 
surface  in  2-D.  A  correction  step  will  probably 
be  necessary  for  curved  boundaries  since  the  con¬ 
dition  staled  is  first-order. 


grid  method,  there  is  still  room  for  improvement. 
Convergence  required  between  100  and  S00  time 
steps,  pointing  out  the  severe  stability  restric¬ 
tions  for  explicit  methods  applied  to  viscous 
problems . 

As  a  final  point,  a  brief  experiment  with  the 

method  of  Leonard  to  difference  the  u  term  was 

x 

conducted,  yielding  excellent  results  (Table  5g,h) 
which  seem  to  verify  the  claimed  third-order  spa¬ 
tial  accuracy  of  the  method.  Otherwise,  central 
differences  on  the  convective  terms  were  used  in 
both  the  1-D  and  the  2-D  viscous  cases. 


l)i  f  f  e  reuc  i  ng 

As  pointed  out  above,  exclusive  use  of  central 
differencing  in  solving  the  grid  equations  may  re¬ 
sult  in  numerical  difficulties.  A  novel  idea  to 
replace  the  usual  "arithmetic  mean  difference," 

Xj.  =  <V\  +  Axl/2A£  with  a  "geometric  mean  differ¬ 
ence,"  -  /  VxAx/AC  was  tested  and  found  to 

i\<  <*  l  I  ini  l  h'Milt  s  in  l-l)  (including  the  1-1) 
wave  pi op.igai  i on  results  presented).  However,  the 
idea  was  disc. tided  for  general  use  because  exten¬ 
sion  to  higher  dimensions  is  difficult.  One  prob¬ 
lem  is  that  Vx  and  Ax  may  be  of  opposite  sign. 
Also,  the  geometric  difference,  operator  is  nonli¬ 
near,  and  the  computation  time  required  to  take 
square  roots  becomes  significant  in  2-D. 


KKSULTS  AM)  DISCISSION 


Two-Dimensional  Inviscid 

Solutions  to  u  ♦  cu  +  du  =0  were  computed 
t  x  y  r 

with  wave  speeds  c  =  d  =  1,  an  initial  condition 
of  a  1  -  1/2  -  0  discontinuity  along  the  main  di¬ 
agonal,  and  a  16x16  grid  on  the  unit  square.  In¬ 
flow  boundaries  were  treated  by  simulating  an  in¬ 
finite  domain  with  a  diagonal  wave,  while  the 
outflow  boundaries  again  required  no  special 
treatment.  Grid  plots,  including  lines  of  con¬ 
stant  u,  are  presented  for  successive  time  steps 
(Figs.  3a-d)  for  a  test  run  with  a  clustering  con¬ 
stant  of  a  =  10.  Good  tracking  of  the  wave  was 
obtained,  and  wave  speed  errors  were  small.  Clus¬ 
tering  at  the  discontinuity  produced  skewed  grid 
cells  that  may  or  may  not  be  desirable  in  more 
general  problems.  Dissipation  in  the  solution  u 
did  occur  to  a  significant  exl. *nt  as  time  pro¬ 
gressed,  causing  the  clustering  of  the  grid  to  di¬ 
minish.  Numerically,  some  problems  were  experi¬ 
enced  when  the  wave  intersected  the  corners, 
forcing  the  use  of  a  small  first  time  step. 


One-li  mimmi*.  i oii.i  1  1  nv  i  sc  id 

Adaptive  grid  solutions  to  u  +  uu  =0  were 

LX 

((imputed  for  the  initial  condition  of  a  1-0  dis¬ 
continuity  on  an  11  point  grid  covering  the  inter¬ 
val  |u,  t  i  .  The  left  boundary  was  fixed  at  u  5  1 , 
while  tin*  right  boundary  required  no  special 
treatment  due  to  the  use  of  upwind  d i f f erenc i ng . 
Results  an*  presented  in  the  form  of  a  grid  time 
lusLniy  plot  (Figs.  2<i,b).  it  can  be  seen  that 
the  grid  tracked  the  wave  quite  well  and  main¬ 
tained  a  reasonable  degree  of  clustering.  Wave 
speed  errors  (compared  to  the  exact  speed  of  1/2) 
were  calculated  and  found  not  to  exceed  2%.  Prob¬ 
lems  with  dissipation  in  u  were  not  experienced, 
partly  si  rue  the  nonlinear  nature  of  Burgers’ 
equation  causes  profiles  having  a  negative  slope 
to  steepen . 


One  -D  uii'Mis  i  on  a  I  V  i  scans 

Steady-state  solutions  to  u  +  uu  =  uu  were 

t  X  XX 

computed  on  a  21  point  grid  for  the  boundary  con¬ 
ditions  u(0)  '  I  and  u(l)  n.  The  initial  dis¬ 
tribution  of  u  was  a  ramp  function  given  by 
u(xj  ~  1 -x  Results  for  viscosity  coefficients  of 
M  =  i.i.nr)  and  u  ~  0  in  are  presented  in  Tables  5 
a-f  and  » ompared  with  the  exact  small  viscosity 
solution,  * . I x  )  -  t anh( ( 1 -x  )/2u) .  Although  the 
adaptive  grid  method  does  provide  better  boundary 
layer  resolution  and  accuracy  than  a  fixed  uniform 


Two-Dimen* tonal  Viscous 

Steady-state  solutions  to  u  +  u  +  u  = 

t  x  y 

p(uxx+u  )  were  computed  on  a  16x16  grid  with  a 

viscosity  coefficient  of  y  =  0.10.  The  boundary 

conditions  were  u(x,0)  =  u,  (x)/u,  (0), 

b  b 

u(0,y)  =  u^Cyl/u^fO),  u(x,l)  =  u(l,y)  =  0,  where 
the  function  is  defined  by  t»b(s)  = 

1  -  exp[ (s-l)/y ] .  The  initial  distribution  of  u 
was  obtained  as  indicated  in  the  preceding 
"Smoothing"  section.  An  exact  t ime- invar iant  so- 
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lution,  u(x,y)  =  ub(x)ub(y ) / ( ub(0) ]  ,  is  presented 

by  Anderson  and  Rai.12  The  computed  results  show 
maximum  deviations  from  this  solution  of  0.019 
(uniform  square  grid),  0.016  (a  =  10),  and  0.014 
(a  =  20).  Resolution  of  the  boundary  layer  was 
judged  to  be  adequate  but  not  ideal.  Convergence 
required  up  to  400  time  steps. 


Cluster ing 

A  major  concern  in  adaptive  grid  work  is  that 
of  obtaining  adequate  resolution  of  the  physical 
domain  with  as  few  total  grid  points  as  possible. 
The  clustering  function  used  in  this  study  is  de¬ 
rived  from  an  equidistr ibut ion  law  based  on  the 
solution  gradient  ,  and  should  provide  reason- 


ible  control  over  grid  point  locations.  Nonethe¬ 
less,  it  has  already  been  demonstrated  that  dis- 
: ret  iz.it  ion  effects  and  smoothing  play  a  crucial 
role  in  determining  the.  grid.  Also,  results  from 
:he  boundary  layer  cases  show  that  achieving  suf¬ 
ficient  resolution  of  viscous  flow  profiles  is  not 
;asy.  For  these  viscous  problems,  it  is  likely 
that  the  solution  gradient  in  physical  space, 

i  -  u,/x,,  should  he  used  to  drive  the  grid  rath- 
x  £  £ 


6T  than  u. 


Such  a  modification  increases  the 


complexity  ot  the  functions  P,  Q  and  their  deriva- 


l  ives  V 


and  Q 


any  event,  it  seems  that  ex- 


T  I 

per  inient  i  i.g  with  alternate  grid  driving  functions 
and  grid  generation  concepts  is  desirable. 


Comput.it  ional  Cons  iderat  ions 

In  one  dimension  the  present  method  is  effi¬ 
cient  in  terms  of  both  computation  time  and  stor¬ 
age,  even  for  large  grids.  Also,  direct  solution 
of  the  simplified  grid  and  grid  speed  equations  is 
poss  i hie  in  1  -I) . 

< in  the  other  h  ind,  stepping  up  to  two  dimen¬ 
sions  greatly  increases  Lime  and  storage  require¬ 
ments.  Direct  solutions  are  no  longer  practical, 
so  iterative  methods  must  be  employed.  The  number 
of  iterative  cycles  required  to  converge  the  grid 
speed  equ.it  ion  may  range  from  only  one,  if  the 
flow  solution  is  not  changing  much  from  step  to 
step,  up  to  30  or  40  (on  a  16xlb  grid)  starting 
from  x^  =  =  0  at  all  points.  Over  40  2-D  ar¬ 

rays  need  to  be  stored  in  order  to  avoid  recalcu¬ 
lating  various  (plant i ties.  Since  explicit  inte¬ 
gration  in  2-D  is  severely  limited  by  small  time 
steps,  it  appears  certain  that  implicit  schemes 
will  have  to  be  developed  for  adaptive  grid  work, 
as  they  already  have  been  for  fixed-grid  PNS  cal- 
culat ions . 


Then*  appears  to  be  little  theoretical  diffi¬ 
culty  in  extending  the  present  adaptive  grid  for¬ 
mulation  to  a  third  dimension,  but  this  might  only 
bo  pi  i'  !  i i  a  1  now  tor  moderate  grid  size  problems 
run  on  a  machine  in  the  supercomputer  class.  Al¬ 
though  adaptive  grids  will  always  he  more  compli¬ 
cated  to  use  than  fixed  grids,  the  extra  work  is 
often  justified  by  improved  solution  accuracy  for 
a  given  number  of  grid  points  ( hence  perhaps  an 
overall  savings  to  achieve  the  same  accuracy  lev¬ 
el). 


As  previously  mentioned,  it  appears  that  the 
present  method  lias  advanced  to  the  point  where  ap¬ 
plication  to  solving  a  set  of  vector  equations  is 
now  possible.  For  invisc id  problems,  an  exp  licit 
pred ictoi /corrector  method*  or  a  split  flux  dif¬ 
ferencing  method  might  !»•»  .».(  t  ess  t  u  1  1  y  employed  to 
i  nt  eg  i  it  e  the  governing  equal  mils.  In  the  viscous 
.  >  I  iss  <»t  non- 1  ter. it  i  ve  ,  implicit,  approxi¬ 
mate  factorization  methods  already  exist  for  solv¬ 
ing  tin’  !’NS  equations  in  generalized  coordinates1 
(mu'  F.i-ipter  M)  However,  interaction  between  the 
grid  and  t low  dynamics  will  force  a  new  look  at 
these  implicit,  schemes  for  adaptive  grid  calcula- 
t  ions 


FUTURE  WORK 

Future  work  should  center  around  the  following 
subjects  already  discussed:  (a)  implicit  integra¬ 
tion  methods;  (b)  curved  or  time-varying  bounda¬ 
ries;  (c)  application  to  the  Euler,  Parabolized 
Navier-Stokes ,  or  full  Navier-Stokes  equations; 

(d)  further  investigation  of  clustering  functions. 
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Tiir  code  employed  herein  has  been  extensively 
ipplifd  t o  ditto  rent  tlow  problems,  including  two- 
i  imi-ns  um.j  1  shock  wive  boundary  layer  interaction 
ner  a  fl.it  plate*,  transonic  turbulent  afterbody 
How1,  turbulent  and  inviscid  transonic  flow  over 
lirtoils11,  and  ot hers  1 ’  1  s .  Turbulence  models  in 
the  examples  above  are  algebraic,  since  completely 
satisfactory  mu  1 1 i -equation  models  for  flows  with 
Jajg"  separated  regions  are  not  yet  available13. 
Du-  v.u  i'-ty  of  tlow  problems  solved  with  this  code 
prove  its  reliability,  making  >t  a  viable  scheme 
to  use  with  this  adaptive  grid  routine. 

INITIAL  Gklll  GENERATION 

Si  me  an  equ  id  i  st  r  i  but  ion  adaptive  grid  rou¬ 
tine  is  not  a  grid  gone rat  ion  routine  as  well,  an 
initial  starting  grid  must  be  created  by  some  oth¬ 
er  means  before  a  solution  is  run.  With  adaption 
available  m  only  one  coordinate  in  the  proposed 
scheme,  the  optimal  starting  mesh  for  this  routine 
would  be  one  which  has  a  sufficient  point  distri¬ 
bution  in  the  remaining  computational  direction. 

Thomas  and  Middlecofi1*  have  introduced  a 
method  of  grid  generation  based  on  earlier  tech¬ 
nique1,  which  allows  for  a  priori  grid  point  clus¬ 
tering  in  at  1  i’js L  one  <  ompui at  iona  1  direction. 

The  technique  follows  the  well-known  method  of 
Thompson,  Thames  and  Mast  in17,  whereby  an  elliptic 
system  nl  two  equal  ions  of  the  form 


r  t  f  =  PK,n) 
’xx  yy 


+  n  =  Q(C,n) 

xx  yy 


is  solved.  By  interchanging  the  roles  ot  the  de¬ 
pendent  (£,n)  and  independent  lx.y)  coordinates  in 
this  equation,  a  quas  i  I  ine.ir  elliptic  system  of 
equations  is  obtained,  which  is  then  solved  by  fi¬ 
nite  differences.  Thomas  and  Middlecofi  have  ob¬ 
tained  analytical  expressions  foi  the  weighting 
functions  P  and  Q  which  will  cluster  points  on  the 
intei  mr  of  the  grid  to  the  same  degree  as  the 
spe<  it  n*d  point  d  i  s  t  r  i  but  ion  on  the  grid  hound- 


TL**  test  grid  described  earlier  (figure  2)  as 
well  as  all  subsequent  grids  were  generated  by  a 
solver  based  on  this  let hn ique .  ]»  Figure  2,  rel¬ 
atively  high  values  lor  were  specified  on  the 

grid’s  inner  boundary  at  the  leading  edge 

(£*■01,)  and  at  both  sides  of  the  trailing  edge 
and  £=hl).  The  influence  from  the  inner 
boundary  is  seen  in  the  point  distributions  for 
n=c  oust  ant  curves  further  from  the  body.  More  im¬ 
port  nil,  though,  is  the  point  distribution  in  the 
n  dii'-ct  itai,  which  does  not  change  with  adaption. 
The  high  grid  clustering  in  the  i;  direction  near 
the  hotly  seen  in  Figure  2  is  necessary  to  resolve 
the  large  solution  gradients  existing  near  the 
body.  With  l  sufficient  \  l  us  ter  mg  in  the 
n -d i roc t  i on ,  the  adaptive  grid  solver  is  then  used 
to  id  lpt  ivelv  cluster  points  in  the  £ -d i rect i on . 

M  MKR  i  ( ]  A I .  KF.Sl’I.TS 

NACA<H»12  Airfoil  Results 

The  majority  of  the  numerical  results  of  this 
sturh  worn  obtained  for  transonic  flow  past  an 
.\A('AMi'l2  airfoil  The  initial  grid  employed  was 
again  the  Otype  2 -D  test  grid  of  Figure  2,  and 
on  this  gi  id.  t  h'*  Euler  equations  were  solved  us¬ 


ing  Tassa’s15  Nav ier -Stokes  code.  The  Euler  equa¬ 
tions  were  selected  rather  than  the  full  Navier- 
Stokes  equations  to  reduce  the  starting  grid 
dimensional  requirements,  thus  permitting  the 
adaption  routine  to  be  tested  on  a  less  complicat¬ 
ed  mesh.  With  the  flow  conditions  selected,  M  = 

0.75  and  a  =  2  degrees,  a  shock  wave  forms  near 
the  upper-surface  midchord  of  the  airfoil. 

A  steady-state  solution  was  run  using  a  vari¬ 
able  time  step  integration  procedure  on  tee  f i  -.ed 
grid.  The  numerical  solution  was  found  to  on- 
verge  after  only  a  few  hundred  time  steps,  result¬ 
ing  in  the  density  contour  field  of  Figure  3a, 
which  gives  no  indication  of  a  shock  wave,  due  to 
the  sparsity  of  grid  points  (Ax-0.6)  in  the  antic¬ 
ipated  shock  region. 

A  new  solution  was  then  formed,  starting  the 
flow  impulsively  from  free  stream  conditions  on 
the  initial  grid,  and  then  passing  control  to  the 
adaptive  grid  routine  after  every  20  time  integra- 
t ion  steps,  again  using  a  variable  time  integra¬ 
tion  procedure.  For  this  case,  minimum  grid  spac¬ 
ing  along  each  n=constant  curve  was  set  at  0.003 

chords  (As  .  =0.003;,  the  initial  weighting  con- 
nun 

stant  A  (equation  (10))  was  set  at  10,  and  both 
functions  f  and  were  turned  off,  allowing  for 

adaption  based  purely  on  density  gradients.  Fig¬ 
ures  8a  and  9a  depict  the  converged  grid  and  den¬ 
sity  field  for  this  test  case.  On  the  newly 
adapted  grid,  the  numerical  density  field  reveals 
a  shock  wave.  The  point  density  in  the  adapted 
grid  has  increased  near  the  shock  region,  but 
still  is  not  exceptionally  high.  As  a  result,  the 
shock  remains  somewhat  smeared,  albeit  over  a 
smaller  region  than  with  the  initial  grid. 

To  investigate  the  convergence  of  this  adap¬ 
tive  grid  algorithm,  it  was  interesting  to  follow, 
among  other  variables,  the  current  value  of  the 
constant  A,  the  actual  minimum  grid  spacing  in  the 
adaptive  coordinate  direction,  and  the  maximum 
distance  any  point  moved  between  adaption  sweeps, 
^Xmax  ^he  adaptive  grid  was  considered  converged 

when  the  first  two  of  these  three  parameters  ap¬ 
proached  a  constant  value,  and  the  last  parameter 
went  to  zero.  The  convergence  history  of  the 
adaptive  grid  test  case  above  is  shown  in  Figure 
7.  These  particular  values  of  the  adaptive  vari¬ 
ables  (A,  B ,  As  .  ,  etc . )  are  seen  to  provide  an 

adaptive  grid  which  converges  quickly  and  smooth¬ 
ly.  Note  that  the  converged  minimum  grid  spacing 

(s(?+l)-8(C))  is  equal  to  0.003,  precisely  the 

value  of  As  chosen  at  the  onset  of  the  problem, 
min  r 

Note  also  that  Ax  reaches  a  peak  value  after 
max 

three  adaptive  sweeps  (60  integration  steps),  and 
not  immediately,  as  one  might  expect.  This  is  due 
to  the  grid  relaxation  technique  defined  by  equa¬ 
tion  (18),  which  prohibits  large  grid  point  move¬ 
ment  in  the  early  numerical  development  of  the  so¬ 
lution  . 

In  order  to  increase  the  shock  region  resolu¬ 
tion  even  further,  it  was  logical  to  reduce  the 

As  .  selected  in  the  adaption  routine,  thus  in- 
min 

creasing  the  final  A  weighting  constant  and  like¬ 
wise  increasing  point  density  in  the  areas  of  high 
fluid  density  gradients.  Figures  8b  and  9b  show 
the  converged  adapted  grid  and  density  contour 
distribution  corresponding  to  a  As^n  equal  to 
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iapt ivo  Algorithm 

ic  numerical  algorithm  used  to  obtain  the  adap- 
ive  grid  point  distribution  governed  by  equation 
9b)  is  reasonably  s t ra lght forward ,  essentially 
ansisting  of  two  rather  large  programming  loops, 
ne  within  another.  The  outer  loop  controls  which 
=  constant  surface  is  being  adapted,  increasing 

roe  n=l  to  The  inner  loop  is  iterative, 

max 

nd  has  a  convergence  criterion  which  must  be  sat* 
sfied  for  the  current  I)  =  constant  curve  before 
he  n  index  is  incremented  on  the  outer  loop.  The 
nner  convergence  loop  works  as  follows:  Equation 
9b)  is  integrated  numerically  to  obtain  a  new 
omi  distribution  £(s).  With  this  vector  and 
i  t,h  p ,  s  ) ,  xlsj  and  y(s),  an  interpolation  scheme 
s  applied  to  calculate  the  newly  adapted  vectors 
namely  s^U).  Pncw<0.  *new«>  and  ' 

he  new  weighting  function  K(s)  (equation  (10))  is 
alculated  from  the  new  vectors,  and  equation  (9b) 
s  oni  n  again  integrated.  This  process  is  contin- 
led  until  the  L  norm  of  the  As  vector,  defined  as 


-2  -  hx  (snew(?)  '  W5))  (l6) 

falls  i>e  low  a  specified  tolerance.  When  the  con* 
/ergeie  e  criterion  is  met,  control  is  passed  to 
the  outer  loop,  ami  the  convergence  loop  is  ap¬ 
plied  to  the  next  n  =  constant  curve. 

At  this  point,  several  comments  about  the  al¬ 
gorithm  are  in  order. 


1.  Due  to  inherent  truncation  errors,  finite  dif¬ 
ference  approximations  of  the  streamwise  density 
gradients  needed  in  function  f^  are  not  suffi¬ 
ciently  smooth.  To  eliminate  this  problem,  it  is 
necessary  to  apply  several  sweeps  of  explicit  sec¬ 
ond  order  smoothing  of  the  form 


=  [1  +  ALP  VHPij 


to  tlie  local  density  field  before  the  grid  is 
adapted,  where  0<ALPS0.25  from  stability  consider- 
at  ions . 

2.  By  adjusting  the  value  of  the  constant  A  in 
equation  (10)  after  each  call  to  the  adaption  rou¬ 
tine,  it  is  possible  to  specify  a  time-asymptotic 
minimum  distance  between  adjacent  grid  points  on  a 

given  s-curve,  called  As  This  constraint  is 

*  m  in 

needed  to  keep  points  from  clustering  too  closely 
across  a  shock,  which  causes  numerical  problems, 
as  explained  later. 

3.  In  many  cases,  particularly  in  the  early  de¬ 
velopment  of  a  flowfield  from  an  impulsive  start, 
numerical  difficulties  may  arise  if  the  grid 
points  are  moved  too  drastically  in  one  adaptive 
sweep  It  is  useful  to  under-relax  the  calculated 
point  d i si r i hut  ion  s(£)  according  to 

S<5)  =  Sol(J(C)h  RELMIN  {  Snew(U  -  SoXd<^} 


where  s(£)  »•  the  final  arc-length  function, 
sold(£)  is  the  current  s  function,  ami  snew^^  18 


the  final  s  function  obtained  from  the  inner  con¬ 
vergence  loop  for  a  given  ti  =  constant  line.  By 
choosing  a  small  value  for  RELMIN  initially,  and 
by  increasing  it  gradually  to  1  after  several 
adaptive  sweeps,  changes  in  the  grid  point  between 
adaptions  are  sufficiently  small  to  prevent  solu¬ 
tion  instabilities. 

4.  (t  is  important  also  to  mention  that  in  the 
present  formulation  of  the  adaptive  algorithm, 
temporal  metric  terms  (grid  speeds)  in  the  trans¬ 
formed  equations  governing  the  flow  are  set  equal 
to  zero.  As  a  result,  once  the  grid  is  updated 
through  adaption,  the  corresponding  solution  vec¬ 
tor  (  e.g,  (p,u,v,e)  )  must  also  be  updated.  This 
is  done  through  interpolation  from  the  current 
grid  and  solution  vector.  Since  the  grid  speeds 
are  neglected,  it  is  difficult  to  determine  if  a 
time  accurate  solution  can  be  obtained  with  this 
adaptive  scheme.  In  light  of  this,  only  steady- 
state  solutions  are  examined  in  this  paper. 


This  algorithm  is  designed  to  supplement  an 
existing  aerodynamic  solver,  ideally  linked  to  the 
main  program  as  a  single  subroutine.  Alterations 
needed  to  implement  the  adaptive  routine  affected 
only  two  per  cent  of  the  total  programming  lines 
in  the  Navier-Stokes/Euler  code  described  next, 
and  it  is  not  anticipated  to  be  much  higher  for 
most  other  codes. 

AERODYNAMIC  SOLVER 

All  numerical  flowfield  results  in  this  study 
were  obtained  from  a  finite  difference  code  devel¬ 
oped  by  Tassa* ,  which  solves  the  unsteady 
2-dimensional  Reynolds  averaged  Navier-Stokes 
equations  written  in  conservation  form  on  a  gener¬ 
al  non-orthogonal  curvilinear  coordinate  system1". 
Flow  variables  and  physical  directions  are  non-di¬ 
mensional  ized  so  that  the  four  governing  P.D.E.s 
and  the  equation  of  state  P  =  pRT,  included  for 
closure  of  the  system,  become  normalized.  This 
allows  the  characteristic  parameters  of  the  flow, 
such  as  Reynolds  number,  to  be  varied  independent¬ 
ly- 

The  resulting  parabolic  system  of  equations  is 
solved  numerically  through  a  modified  form  of  the 
Bri ley-McDonald  Alternating  Direction  Implicit 
scheme11.  Whereas  the  Briley  and  McDonald  dual 
time  level  scheme  represents  all  but  the  energy 
equation  in  conservation  form,  Tassa's  modified 
three  time  level  scheme  writes  even  the  energy 
equation  in  conservation  form.  Non-linear  terms 
in  these  equations  are  linearized  by  using  Taylor 
series  expansions  at  the  known  lime  level.  By 
representing  the  dependent  variables  p,  u,  v  and  e 
as  the  sum  of  values  at  the  known  time  level  and 
an  incremental  value,  a  linear  matrix  equation  is 
obtained  in  terms  of  the  unknown  incremental  val¬ 
ues.  The  Douglas -Gunn1 2  procedure  for  generating 
ADI  schemes  is  then  applied  to  the  new  system  of 
equations,  splitting  the  matrix  equation  into  a 
system  of  two  one-dimensional  matrix  operators. 
After  discretizing  the  spatial  operators  using 
second-order  formulas,  the  incremental  solution 
vector  is  found  by  block  elimination  techniques. 
Tassa  and  Schuster13  have  found  it  necessary  to 
add  artificial  dissipation  near  regions  of  severe 
pressure  gradients  such  as  shocks,  to  suppress 
high  frequency  components.  In  addition,  fourth- 
order  dissipation  is  added  to  the  dependent  vari¬ 
ables  in  the  Euler  equations  in  order  to  reduce 
the  overshoot  of  pressure  across  the  shock. 
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The  density  gradient  is  selected  with  the  case  of 
transonic  flow  in  mind.  Across  a  shock,  density 
changes  rapidly  iti  a  physical  sense,  and  for  the 
Kuler  equations,  discont muons ly  in  a  mathemat ica 1 
sense.  If  the  sLreamwi.se  direction  of  flow  and 
the  -.-curves  (n  =  constant  loci)  are  nearly 
aligned  in  space,  then  3p/ds  will  be  relatively 
high  in  the  shock  region,  forcing  the  grid  points 
to  move  towards  the  shock  location.  For  llio  test 
solnlmn  previously  mentioned,  this  is  indeed  the 
case,  is  evidenced  by  the  adaptive  grid  of  Figure 
3h.  Note  that  along  each  n  -  constant  line,  the 
areas  of  higliei  grid  point,  clustering  are  at  the 
airfoil  leading  edge  and  in  the  shock  region. 

This  is  as  expected,  since,  it  is  in  exactly  these 
regions  where  a  large  st Teamwise  density  gradient 
is  observed  (Figure  3a) 


regard  to  the  points  outside  the  domain,  there  is 
no  guarantee  that  grid  spacing  is  continuous 
across  the  boundary.  Figure  4  illustrates  this 
fact.  Along  the  £  ^  and  £2  1 mes ,  where  grid  spac¬ 
ing  in  the  £  direction  varies  rapidly,  the  accura¬ 
cy  of  the  finite  difference  equations  may  be  inad¬ 
equate.  For  example,  at  any  point  along  the  arc, 
a  finite  difference  equation  with  second-order 
truncation  error  on  a  uniformly  spaced  grid  will 
decay  to  first  order  accuracy  on  a  non-uniform 
grid  whenever  the  grid  spacing  on  both  sides  of 
the  point  differ  greatly-  that  is,  when 
I  As+  |«*>r  I  As+  I  >>1  . 


To  remedy  this  spacing  problem,  a  third  weighting 
function  is  introduced,  defined  as 


f3(s)=  D  Ss(s=smax)e 


-*(r— ) 

\  max/ 


In  oider  to  preserve  the  shape  of  the  s-curves 
while  us ing  adaption,  it  is  impoitant  to  retain  at 
least  a  minimum  number  of  grid  points  in  the  re¬ 
gions  of  high  arc:  curvature,  particularly  along 
the  boundaries  of  the  physical  domain,  where  al¬ 
tering  the  boundary  will  likewise*  alter  the  prob¬ 
lem  und»*r  consideration.  For  this  reason,  f^fs) 

is  del  i  lied  as 

Vu  ‘  Vu 

f .  (s)  =  B  -  =  B  |  K  |  (12) 

<x2+y2,3'2 


'  s  -s 
max 


+  (s=0) e 

s 


D  >  0 ,  g=S0 


where  the  asterisk  refers  to  the  value  of  eval¬ 
uated  outside  of  the  adaptive  boundary.  This 
function  is  chosen  since,  for  a  largo  enough  value 
of  g  (-SO), 


f3(W  5  Df^(smax)1'  f3(0)  * 


dK*(0)] 


where  K  is  i  In*  matliemat  i  ca  I  definition  of  the  cur¬ 
vature  of  s  (parameterized  by  £).  and  B  is  again  a 
positive*  constant,  used  to  control  the  degree  of 
clustering  based  on  grid  curvature.  Observe  that 
if  the  spatial  density  gradients  are  negligible 
and  i  I  both  fun*  t  ions  I  ,  (s)  and  l.^ls)  are  set 

equal  to  zero,  then  tin*  i ight  hand  side  of  equa¬ 
tion  1  Hi )  approaches  unity  With  this  weighting 
function,  gr id  points  are  spaced  at  equal  incre¬ 
ments  .it  tln>  s -d  i  x  <-ct  ion  •iong  each  n  =  constant 
surface,  as  shown  m  figiM*  4.  Compared  with  the 
test  grid  of  Figure*  tin-  grid  in  Figure  4  has 
pooi  resolution  at  the  1  •  *  •  *  *  1 1 1 1  g  edge,  and  hence, 
the  tine  airfoil  shape  is  i,-»i  sufficiently  de¬ 
fine  i  lit  t  t  lie  ntK)i  ,  ',ini  I-  the  sh  ipc*  of  each  n  = 
‘on si  ml  <  urve  r,  ilgebi  ii  aliy  defined  by  curves 
til  through  l he  gi id  points,  few  points  in  regions 
of  high  an.  curvature  might  eventually  alter  the 
shape  of  i  fur  airfoil  j  f  t  *•  i  several  adaptive 
swei-ps .  Fortunately,  by  c  lusteiing  points  in 
these  high  curvature  regions,  this  problem  can  be 
a  vo  i  led.  Figure  ’>  depicts  the  same  test  grid  with 
points  redistributed  awonliug  lo  function  f^  for 

1  notniNu.  value  ot  the  constant  H.  (ju  this  grid, 
t  he  ‘.I,  ipe  of  the  t)  -  constant  curves  near  the 
leading  edge  are  preserved, 

Another  noticeable  dilteieiue  between  Figures 

2  and  w  is  in  the  grid  spa:  ing  tieir  the  boundaries 
ot  i  lie  ad<ipLive  domain.  The  elliptic,  properties 
of  t  lie  ojiginal  grid  generator  insure  tfiat  grid 
spa<  ing  c fianges  -unooLhlv  t  h rough* ml  the  mesh 
However,  since  the  adaptive  grid  algorithm  redis- 
tiihutes  points  within  the  adaptive  domain  with  no 


f  3  (s>  *  0  for  0«  S<<  smax 


(14a. b) 


Differentiating  equation  (9b)  whith  respect  to  s 
yields  /- 

Is  max 

F  (s)  =  5-  JO  F(s)ds  =  5s<ci>  05) 


( 


Provided  that  both  grid  curvature  and  density  gra¬ 
dients  are  negligible  at  the  adaptive  domain 
boundaries,  it  can  be  shown  that  by  choosing  the 
constant  D  =  C^/A,  grid  spacing  (£s>  will  be  con¬ 
tinuous  across  the*  boundary  of  I)  ,  As  a 

adapt ive 

result,  the  order  of  the  solution  truncation  error 
near  the  boundary  should  not  decrease  (still 
2 

o(Ax  )  )  appreciably. 

Function  f,^(s)  is  used  to  create  the  adaptive 

grids  in  Figures  6a  and  6b.  The  first  of  these 
grids  has  a  weighting  function  equal  to  f^(s) 

alone,  producing  a  grid  with  equal  spacing  along 
each  s-r.urve  everywhere  except  near  the  boundary 

of  D  ,  Both  functions  f_  and  f  are  used 

adapt ive  2  3 

to  create  Figure  6b,  which  has  clustering  near  the 
leading  edge  as  well  as  equal  grid  spacing  through 
the  boundary.  Indeed,  it  is  a  combination  of  all 
three  weighting  functions  that  will  produce  the 
most  desirable  adaptive  grid. 


Consider  first,  a  positive  weighting  function 
w=w(s),  associated  with  some  partial  differential 
equation,  chosen  so  that  w  increases  as  the  grid 
point  density  (  £s  )  needed  to  approximate  the  so¬ 
lution  to  the  partial  differential  equation  to 
some  fixed  error  also  increases.  Saltzman*  has 
shown,  through  a  variational  approach,  that  by 
minimizing  the  integral 

s 

/max 

W(s)  d£  (5) 

0 

the  error  due  to  solution  approximation  is  also 
minimized.  The  Eu ler-Lagrange  equation  corre¬ 
sponding  to  this  integral  is  now 


For  a  general  weighting  function  F(s),  neither 
of  equations  (0)  can  be  nnmerical ly  integrated  di¬ 
rectly.  The  integrands  on  the  right  hand  side  of 
each  equation  are  functions  of  the  grid  point  dis¬ 
tributions  s(£)  and  £(s),  which  appear  also  on  the 
left  hand  side  of  each  equation.  To  solve  the 
equations,  then,  an  indirect  method,  such  as  an 
iterative  updating  procedure,  must  be  employed. 
Although  equation  (9a)  may  appear  to  be  a  better 
choice  than  (9b)  for  numerical  integration,  on  the 
basis  of  iterative  convergence  speeds  this  was  not 
the  case.  Rather,  the  number  of  iterations  needed 
for  convergence  of  (9b)  was  as  much  as  an  order  of 
magnitude  less  than  the  number  needed  for  conver¬ 
gence  of  (9a),  and  because  of  this,  equation  (9b) 
was  used  exclusively  in  this  work.  This  differ¬ 
ence  in  convergence  speeds  was  due  in  part  to  the 
form  of  the  weighting  function  F(s)  used,  de¬ 
scribed  below. 


/  _3 _ 3  _3_  \ 

V  3 f  3s  Hs  /<W  S£  *  =  0 

whuli  induces  to  tbo  condition  that 


(6) 

Weighting  Function 

The  weighting  function  chosen  for  this  study 
is  of  the  form 


W(f.l 


constant , 


(7a) 


F  (s)  =  1  +  AJfjIs)  +  f2 


(s)  +  f3<s)) 


(10) 


or,  with  s  as  thft  independent  variable, 

W(s)*5  /  £  =  constant  (7b) 


since  s  is  assumed  to  be  a  function  of  £  only. 
Replacing  /w(£)  wiih  F(£)  for  the  sake  of  conven¬ 
ience,  and  assuming  that  the  A£  between  adjacent 
grid  points  along  each  s-curve  is  equal  to  one, 
equation  (7a)  can  be  approximated  by  forward-dif¬ 
ferences  as  F(£ )As=constant ,  where  As 
=s(5+l )-s(C) .  This  states  that  the  product  of  the 
weighting  function  and  the  grid  spacing  is  equally 
distributed  along  the  s-curve.  Anderson  et  al.' 
have  for  this  reason  named  adaptive  grid  methods 
governed  by  equation  (7)  equ id i s t r ibut ion  schemes. 
Now,  the  boundary  conditions  for  (7a)  and  (7b) 
along  each  s-curve  are 


S(F.j)  =  0 

s  ( c )  =  s 
2  max 

and 


(8a) 


U0)  = 

Usmax>  =  C2 


(8b) 


The*  i-orresponding  solutions  to  equations  (7a)  and 
(7b)  nro  then  found  to  be 


and 


1/F<£)  d£ 


1/F(£)  d£ 


(s) 

max 

(s) 


ds 


ds 


(9a) 


(9b) 


respect i ve 1 y . 


where  A  is  a  positive  constant  and  fj(s),  f2(s), 
and  fj(s)  are  each  non-negative  functions.  In¬ 
cluding  the  constant,  1,  in  F(s)  allows  A  to  con¬ 
trol  the  degree  of  grid  clustering,  and  insures 
that  F(s)  will  not  approach  zero  (£^0),  which  is 

not  feasible. 

To  illustrate  the  utility  of  each  of  these 
three  functions  fj,  f^  and  f^,  a  test  grid  was 

generated,  the  inne.r  detail  of  which  is  presented 
in  Figure  2.  The  inner  boundary  of  this  C-type 
grid  (n  =  1  surface)  is  an  NACA0012  airfoil,  with 
61  points  wrapped  counter-clockwise  around  the 
surface,  and  92  more  points,  ranging  from  )££S21 
and  81££S101,  distributed  downstream  of  the  trail¬ 
ing  edge.  The  adaptive  domain,  Da<japtjve> 

bounded  by  two  inconstant  lines  emanating  from  the 
trailing  edge  (£ j=2 1 ,£2=8 1 ) .  The  adaptive  coordi¬ 
nate  chosen  is  £,  meaning  that  s(£)  functions  will 
be  redistributed  along  each  of  the  21  nearly  con¬ 
centric  inconstant  curves.  The  dimensions  of  the 
adaptive  and  total  domains  are  then  61  x  21  and 
101  x  21,  respectively.  On  this  grid,  a  converged 
2-D  conservative  variable  Euler  equation  solution 
was  generated  for  a  Mach  number  M  =0.75  and  angle 

of  attack  a  =2 . 0  degrees.  These  flow  conditions 
are  known  to  produce  a  shock  just  upstream  of  the 
upper  surface  midchord  oil  a  grid  with  sufficient 
point  density  in  the  shock  region.  On  this  test 
grid,  however,  the  shock  region  point  density  is 
sparse,  and  the  shock  is  smeared  across  several 
grid  points,  as  shown  in  Figure  3a,  which  pictures 
curves  of  constant  fluid  density.  The  methods 
used  to  generate  both  the  initial  grid  and  the 
initial  solution  are  described  later  in  this  pa¬ 
per. 

The  most  important  term  in  the  weighting  func¬ 
tion  is  fj(s),  defined  as  the  first  partial  deriv¬ 
ative  of  fluid  density  with  respect  to  arclength 
s.  Numerically,  this  derivative  is  easily  calcu¬ 
lated  by  noting  that 


3 


minimizing  the  integral,  additional  P.D.E.s  must 
be  solved.  Nevertheless,  the  fact  that  these 
techniques  are  based  on  a  firm  mathematical  foun¬ 
dation  should  eventually  make  them  more  popular 
than  ad  hoc  procedures. 

Ono  particular  type  of  variational  adaptive 
grid,  previously  employed  by  both  Dwyer*  and  Gnof- 
fo\  is  particularly  attractive  if  adaption  is 
needed  in  only  one  computational  coordinate.  This 
technique  is  referred  to  as  an  equ idist r ibut ion 
adaptive  grid  scheme,  and  is  the  type  of  grid 
scheme  extensively  studied  in  this  work.  The 
scheme,  which  is  formulated  and  explained  in  some 
detail  in  the  next  section,  is  an  ideal  technique 
for  use  with  transonic  airfoil  problems,  where 
grid  point  adaption  is  usually  only  needed  in  the 
streamwise  direction  of  flow. 


ADAPTIVE  GRID  SCHEME 


Mathematical  Formulation 


When  f in i te-di f terence  techniques  arc  used  to 
obtain  the  solutions  to  the  partial  differential 
equations  governing  fluid  flows,  a  finite  number 
of  points  in  physical  space  must  first  be  select¬ 
ed.  These  points  comprise  the  grid,  or  mesh,  at 
which  the  solution  to  the  discretized  versions  of 
the  P.D.E.s  (the  finite  difference  equations)  are 
to  be  calculated.  In  two  dimensions,  the  physical 
location  of  each  grid  point  can  be  defined  by  its 
two  Cartesian  coordinates,  (x,y).  Mesh  points  can 
aiternat ively  be  defined  by  two  coordinates  4  and 
n,  chosen  so  that  the  grid  points  in  the  x-y  phys¬ 
ical  plane  become  equally  spaced  and  fixed  in  time 
in  the  4_n  computational  plane.  For  simplicity ,  4 
and  n  are  integer  valued,  usually  set  equal  to  1. 
The  computai iona 1  coordinates  of  each  point  are 
then  associated  with  a  storage  location  in  a  two- 
dimensional  array.  The  representation  in  either 
the  computational  or  physical  plane  uniquely  de¬ 
fines  a  given  mesh  point. 


Now,  the  computat iona 1  plane  and  the  physical 
plane  are  mathematically  related  through  the  vec¬ 
tor-valued  mapping 


M 

x  ( f, ,  n ) 

M 

y  (C,n) 

(la) 


which  is  schematically  represented  in  Figure  1. 
Provided  that  this  transformation  is  both  one-to- 
one  and  onto,  the  inverse  mapping  also  exists,  and 
is  defined  as 


N  , 

5 (x,y) 

rl 

n (x,y) 

(lb) 


In  differential  notation, 
be  written  as 

/  i  r 


this  transformation  can 


<34 


(2a) 


dy 


dn  J 


or  inversely  as 


d£ 

y  -x 

*  n  n 

dx 

.  I 

•  1  1 

(2b) 

dn 

;yt  xt. 

dy 

whom  J  *=  x,y  *x  y.  *  the  Jiicnhinn  of  the  Mapping. 

S  h  1  H 

The  major  advantage  of  employing  coeiputat  iona  I  co- 
ordinates  is  thaL  through  equations  (2a)  and  (2b), 
the  P.D.E.s  governing  the  flow  can  be  transformed 
so  the  independent  variables  are  now  £  and  q . 

This  reduces  the  complexity  of  the  finite-differ¬ 
ence  equations,  since  the  computational  grid  on 
which  they  are  solved  is  both  equally  spaced  and 
non-moving. 


The  complete  grid  region  in  the  physical  plane 
will  be  known  as  the  total  domain,  and  is 


mathematically  defined  as 
/ 

D 

x  -  x((,n) 

1  <  r  <  r 

total « 

y  =  y(C»n) 

1  <  n  <  n 
*  -  1  max 

The  ranges  on  the  computational  coordinates  are 
chosen  for  convenience  only.  Either  past  experi¬ 
ence  or  intuition  may  determine  that  inadequacies 
in  grid  point  spacing  prevent  the  discretized 
P.D.E.s  from  approximating  the  governing  equations 
of  motion  to  a  desired  degree.  In  general,  only  a 
subregion  of  the  total  domain  will  have  an  inade¬ 
quate  grid  point  distribution,  so  only  the  grid 
points  in  this  region,  called  Dadaptive>  and  de" 


fined  as 

Dadaptive 


l  x  *  x (t ,n) 

[  y  -  y(S»n) 

(4) 

will  need  to  be  redistributed.  D  ,  as  the 

adaptive 

name  suggests,  is  then  the  domain  in  which  the 
adaptive  transformation  will  be  applied.  Schemat¬ 
ically,  ®adaptjve  *s  the  shaded  area  of  Figure  1. 

Note  that  since  DadaptlveC  D^j,  the  cases 

where  either  {,“1  or  (  =(  are  certainly  a 1 Iowa- 
x  l  max 

ble . 


Now  let  s({)  be  the  arclength  along  an 

inconstant  surface  in  D  .  such  that 

adaptive 

«({=(,  )“0  and  s  (£«=(_  )=s  _  ,  and  observe  that 
*  4  max 

s(l)<s(5+l)  for  all  .  Calculation  of 

a  new  s(()  function  for  each  inconstant  surface 
will  result  in  a  redistribution  of  the  points 
along  the  s -curve.  Since  the  points  are  free  to 
move  only  along  inconstant  surfaces  and  not  along 
("constant  surfaces,  the  selection  of  the  adaptive 
coordinate  ((  in  this  case)  is  not  a  trivial  mat¬ 
ter,  but  Is  dependent  on  the  type  of  flow  being 
modeled.  Fortunately,  in  many  cases,  particularly 
those  with  well-defined  regions  of  lsrge  solution 
grsdients,  the  correct  adapting  coordinate  is  im¬ 
mediately  apparent.  The  problem  now  becoaes  one 
of  determining  the  new  s(()  function  for  each 
inconstant  curve. 
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ABSTRACT 

A  dynamically  adaptive  grid  scheme  based  on 
equidistr ibut ion  in  one  computational  coordinate 
is  applied  for  the  first  time  to  inviscid  transon¬ 
ic  flow  numerically  solved  on  C-type  airfoil 
grids.  Steady-state  solutions  arc  obtained  for 
NACA0012  and  RAE2822  airfoils  using  both  fixed  and 
solution  adaptive  grids,  and  results  for  both 
grids  are  compared  with  previous  numerical  and  ex¬ 
perimental  data.  The  adaptive  grid  algorithm  is 
seen  to  resolve  details  of  the  flow  field  near  the 
upper-surface  midchord  shock  not  seen  in  the  fixed 
grid  solution,  thus  eliminating  the  need  for  a 
priori  grid  point  clustering  in  the  region  of  the 
anticipated  shock.  In  addition,  problems  inherent 
to  schemes  of  this  type  are  discussed,  and  sugges¬ 
tions  for  further  study  are  also  made. 


INTRODUCTION 

The  generation  of  compnl.it  iona  I  meshes  for  use 
in  so  King  d  i. sc  n- L  i  zed  systems  of  pari  ini  differ¬ 
ential  equations  fP.D.K.s)  is  presently  a  subject 
of  intense  research.  In  most  cases,  a  grid  is 
generated  by  either  algebraic,  complex  variable, 
or  differential  equation  methods  before  any  numer¬ 
ical  solutions  are  calculated,  with  the  resultant 
mesh  being  used  for  all  subsequent  computations. 

A  major  pitfall  of  this  accepted  technique  lies  in 
the  fact  that  the  mesh  point  distribution  often 
proves  to  be  inadequate  for  approximating  the 
problem  under  consideration. 

When  the  P.D.E.s  governing  certain  fluid  flows 
are  discretized  and  solved  on  an  insufficient 
mesh,  it  is  possible  that  certain  high  gradient 
phenomena  of  the  flow  field  will  not  be  captured, 
due  to  the  sparsity  of  grid  points  in  the  high 
gradient  regions.  For  example,  boundary  layer, 
free  shear  layer,  and  captured  shock  regions 
(flows  with  multiple  length  scales)  all  require 
locally  high  grid  resolution  for  the  solution  to 
be  approximated  to  a  given  degree  ot  accuracy. 

This  problem  might  be  alleviated  by  using  a  grid 
with  high  resolution  throughout,  but  the  added 
computer  storage  and  time  demands  make  thia  im¬ 
practical.  Alternatively,  points  could  be  clua- 
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tered  only  in  the  anticipated  regions  of  large 
solution  gradients.  This,  however,  requires  an  a 
priori  knowledge  of  the  flow  field,  which  many 
tines  is  not  available. 

A  better  suggestion  is  to  use  a  dynamically 
adaptive  grid,  i.e.,  one  which  continually  adjusts 
the  point  distribution  within  the  mesh  as  the  flow 
solution  is  advanced  in  tine.  An  ideal  adaptive 
mesh  scheme  would  be  one  which  readjusts  the  grid 
points  so  that  the  local  truncation  errors  of  the 
discretized  equations  of  motion  are  reduced  to  a 
minimum,  constant  value  throughout  the  mesh.  Such 
s  scheme  would  require  analytical  or  approximate 
expressions  for  the  local  truncation  error.  Un* 
fortunately,  except  for  the  most  basic  equations 
of  fluid  motion,  expressions  for  the  truncation 
erroi  of  the  associated  finite  difference  equa¬ 
tions  are  extremely  difficult  or  impossible  to 
calculate.  Therefore,  when  truncation  errors  are 
not  available,  grid  point  adaption  should  be  driv¬ 
en  by  some  other  mathematical  or  physical  relation 
which  will  still  improve  the  overall  accuracy  of 
the  finite  difference  equations  of  motion.  For 
most  adaptive  grid  schemes,  this  is  indeed  the 
case,  although  the  means  used  to  reach  this  end 
differ  greatly. 

Anderson'  has  separated  existing  grid  tech¬ 
niques  into  two  distinct  categories.  In  the  first 
of  these  :wo  categories,  a  mathematical  law  defin¬ 
ing  the  speed  of  the  grid  points  is  postulated, 
and  the  new  grid  point  locations  are  obtained  by 
integrating  in  time.  Although  calculating  grid 
point  locations  from  the  grid  speeds  is  straight¬ 
forward,  formulating  grid  speed  laws  based  on 
sound  physical  reasoning  is  not.  In  fact,  a  fair 
amount  of  Ingenuity  is  often  necessary.  For  exam¬ 
ple,  in  a  technique  developed  by  Rai  and  Ander¬ 
son*,  grid  speeds  are  calculated  from  an  attrac¬ 
tion  model.  In  this  model  every  two  grid  points 
induce  on  each  other  a  small  velocity  which  is  de¬ 
pendent  both  on  the  distance  between  the  two 
points  and  the  other  point's  deviation  from  the 
average  solution  error.  The  grid  speed  at  each 
point  is  then  set  equal  to  the  sum  of  all  of  the 
velocities  induced  by  every  other  point  in  the 
grid.  A  survey  of  adaptive  grid  schemes  based  on 
grid  speed  laws  is  given  in  reference  3. 

In  the  second  of  these  two  categories,  grid 
points  are  moved  to  new  locations  through  speci¬ 
fied  mathematical  mappings,  and  then  the  grid 
speeds  needed  in  the  equations  of  motion  are  cal¬ 
culated  from  backward  differences.  Among  the 
techniques  falling  into  this  class  are  variational 
methods,  developed  extensively  by  Brackbill*  and 
S — 1 1 m an  In  variational  methods,  an  integral 

containing  a  measure  of  some  grid  parameter,  such 
as  orthogonality  or  smoothness,  is  minimized.  One 
disadvantage  of  variational  method*  is  that  by 
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Table  5.  One-dimenslonal  viscous  Burgers'  equation.  Numerical  results  for  boundary  layer. 


(a)  u  *  0.05.  a  =  0.0 

step  =  20 0,  time  *  5.00 


20  O.'/Hu  0.2911  0.0094 

21  1.0000  0.0000  0.0000 

Integrated  error  -  0.00129 


(c)  ■■■  O.O'i,  a  ;’0.0 

step  =  010,  time  ,’./2 


N 

X 

a 

error 

1  / 

0.8815 

O.H44  1 

0.0153 

IK 

0.9186 

0 .  1)8  7 1 

0.0154 

19 

0.9489 

0.48  12 

0.0122 

20 

0.9753 
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0.0000 

0.0000 
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0 
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-  3.7’> 

\ 

,, 

error 

1  5 

0.9109 

0 . 0057 

K. 

0.  /V)(» 

0.83)/ 

0.0074 

1  ) 

0.800)) 

0.  //ill, 

0.0090 

IK 

0.8500 

0 . 1,3  48 

0.0096 

l*» 

O.OMmii) 

0.  70t, 

0.0085 

21) 

0.95(H) 

>t)0 

0.0051 

21 

i  .oooo 

o.ot  too 

o.oooo 

1  nt  eg 

’.•lull  i-rr<> 

r  -  0.002 

96 

(e)  u  -  0.10,  a  «  10.0 

step  -  400,  time  »  3.43 


N 

X 

u 

error 

17 

0.8000 

0.9756 

0.01 16 

18 

0.8500 

0.9286 

0.0234 

19 

0.9000 

0.8000 

0.0184 

20 

0.9500 

0.5000 

0.0379 

21 

1 . 0000 

0.0000 

0 . 0000 

Integrated  error  3  0.00601 

(b)  „•  =  0. 

05,  <1  *  10. 

0 

step  = 

:  300,  time 

1 .  19 

N 

X 

u 

error 

17 

0.8570 

0.9000 

0.0156 

18 

0.9027 

0. 

0.0188 

19 

0.9 

</.  Oh 2 

0.0164 

N 

X 

u 

error 

15 

0.7388 

0.8687 

0.0056 

16 

0.7887 

0.7906 

0.0064 

17 

0.8360 

0.6816 

0.0066 

18 

0.8801 

0.5419 

0.0060 

19 

0.9219 

0.3765 

0.0045 

20 

0.9614 

0.1929 

0.0024 

21 

1.0000 

0.0000 

0.0000 

Integrated  error  -  0.00224 

(f)  V  *  0. 

10,  a  =  20 

.0 

step  « 

500  time 

=  3.53 

N 

X 

u 

error 

15 

0.7624 

0.8354 

0.0054 

16 

0.8098 

0.  7460 

0.0058 

17 

0.8532 

0.6310 

0.0056 

18 

0.893) 

0.49  )6 

0.0048 

19 

0.9302 

0.3389 

0.0035 

20 

0.9656 

0.1724 

0.0018 

21 

1 . 0000 

0.0000 

0.0000 

Integrated  error  =  0.00204 


(g)  Leonard  Ird-order  differencing 
0  =  0.05,  a  =  5.0 
step  -  100,  time  =  3.21 


X 

u 

error 

0.7810 

0.9789 

0.0037 

0.8826 

0.8394 

0.0137 

0.9501 

0.4867 

0.0251 

1 . 0000 

0.0000 

0.0000 

Integrated  error  *  0.00312 


(h)  Leonard  3rd-order  differencing 
U  *  0.05,  a  =  5.0 
step  =  200,  time  =  2.66 


N 

X 

u 

error 

16 

0.7846 

0.9728 

-0.0006 

17 

0.8364 

0.9269 

-0.0001 

18 

0.8858 

0.8177 

0.0025 

19 

0.9292 

0.6157 

0.0061 

20 

0.9663 

0.3306 

0.0054 

21 

1 .0000 

0.0000 

0.0000 

Integrated  error  =  0.00047 


Table  3.  Finite  difference  calculation  of  P^  and  QT 


Note:  Examples  given  for  calculation  of  P  .  Expressions  for  are  slmila 
General  form 


P  -  P(x,y,u) 


P:lj  =  C,P1. /U,k.  l)U.k  ,  +  (;)Pl.J/3Xk.l)Xlkl  +  (3Pl.j/3yk.l,yi 


where  summation  on  the  indices  k,l  is  Implied. 


Specif ic  form 

P(x,y,u)  -  P(u)/Y(x,y) 


PT  -  -pyt/y  +  pt/y 


\  -  3/9t  <xf2  +  *?>  3  2(x5xtc  +  Vxt) 


(Yp  ,/3u  )u  as  above. 


i.J 


i,  j  k,l  i 


k,l 


For  P  *  (2aufuf  f  )/(l  +  au{.)  with  central  differencing, 


(',Pi,j/:Vl.j)  ’  P(2VU1.J  +  Fi.j\tJ) 


OPi.J/aUi.J)  “  '4pUf. 


I.J 


(aPi,j/'<Ui+l,J)  ~  p(2Aui,J  *  Pi.JU5i  ^ 


where  p  *  a/(l  +  au,  ) 


i.J 


k.l 


I  able  4.  Smoothing  operator 


Table  1.  (Cont.) 


1-D  equations 


8.  Governing  equation 


2-D  equations 


(a)  General  conservative  fora 


u.  +  f  -  viscous  terns 
t  x 


ut  +  £tue  +  cxfc 


=  viscous  terms 


(b)  Burgers'  equation 


+  UU  a  UU 
t  X  XX 


UT  +  +  U^)u£  a  JlV  U 


+  f  +  g  -  viscous  terms 
t  x  y 


UT  +  ?tu5  +  ntun  +  Vc  +  nxfH 

+  £  jr  +  n  g  *  viscous  terms 

y  {  y  n 


u^  +  cu  +  du  ■  uV  u 

t  x  y 


UT  +  (St  +  C5x  +  V? 


+  <nt  +  cnx  +  dny)un  »  pv  u 


See  equation  set  (2)  for  V  .  Note  that  VC*?  and  V  r|  »  Q. 


Table  2.  Coefficient  matrices  [A]  and  [B]  in  2-D  grid  speed  equation 


(A)  -  (ar  a2I  [B]  -  [bi.  b2J 


Let  r  =  (x,y)  and  k  -  J(Pr^  +  Qr^) 


a,  -  2(-x  rf  +  xrr  +  y  k) 

1  n  f,n  C  nn  n 

*■  -*-*-*• 
a.  *  2(-y  rf  +  y.r  -  x  k) 

2  n  Cn  C  nn  n 


bi  ‘  2(Vct  ’  VCn  -  yCk) 
b2  ‘  2(VcS  ’  yC"Cn  +  xCk) 


Table  1.  (Cont.) 


1-D  equations 


2-D  equations 


Steady  grid  equation 

(a)  V2£  =  P(C.u) 

(a)  V2e  -  P(S.n.u)  ;  v2n  -  Q(C.n,u) 

(b)  x^r  +  x2P(x,u)  *  0 

(b)  ara  -  2dr^  +  Y?nn  +  J2^  +  <£„> 

where: 

r  -  (x,y)T  6  -  xfxn  +  y^ 


4.  Integral  grid  law 


x 

J  [ 1  +  aw(s) Ids 
- - - r 

,  ’max 

/  [ 1  +  aw(s) ] ds 

0 


5.  Clustering  function 
awf 

I*  <  x .  ii  /  »  - 

X.  ( I  +  .iw) 

where: 
w(u)  ■=  uf 


p(x>y-u)  ■  w T7aw;y 


av2C 

Q(x’y‘u)  "  a(l~\w2) 


where: 

2  2 
wx(u)  »  W2(u)  - 

a,y  defined  in  (3). 


b.  Simplified  grid  equation/clustering  function 

P(x,u)  =  P(u)/x2  P(x,y,u)  -  P(u)/y(x,y) 

x.f  +  xrP(u)  •  0  Q(x,y ,u)  -  Q(u)/a<x,y) 

Substitute  into  (3b). 


Grid  speed  equation 
(,l)  (x|  ).  .  + 

+  x?P  =0 

.  I 

(b)  (xT)fJ.  +  P(xr)f  +  x.PT  -  0 


(a)  a(rT)5f;  -  +  Y(rT>nr) 

+  [A]UT)C  +  [B)UT)n 

+  j2<V5  +  vn>  ■  0 

See  l A]  and  [B]  in  Table  2. 
See  P  ,  Q  in  Table  3. 
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Fig.  4  Two-dimensional  viscous  Burgers'  equation 
converged  grid  plot. 

(a)  a  -  10.0,  steps  -  300,  time  ■  1.90 

(b)  a  -  20.0,  steps  -  400,  time  -  1.98 


Table  1.  Summary  of  equations 


1-1'  equations 


1.  Coordinate  transformation 


=  f.(x,t) 


2~D  equations 


f.  »  C(x.y.t) 
n  -  n(x,y,t) 


pnHHHHI 

MiMieiBMlHH 
HBSBmm— ■■■ 

msaa 

MM 


■gsaa 

BBSS^ 

■■■■■■■■■qSsSSSI 

■■■■■■■■■■■■■Si 


LT-IDTH 

■■■■■■■■■bBSsBS 

■■■■■■■■■■bSSSS 

lilBBSBBiSi"S8S 


Two-dimensional  lnviscid  Burgers  equation 
grid  plot  (a  -■  10.0,  u)  -  1.0,  X  •  20.0). 

(a)  stepB  «  0,  time  •  0.00 

(b)  steps  -  5,  time  ■  0.05 

(c)  steps  -  10,  time  -  0.10 

(d)  steps  -  15,  time  -  0.15 

Note:  Dashed  lines  define  boundary  of  wave 
and  center  of  wave. 


0,0025.  As  expected,  the  shock  region  in  these 
figures  is  better  resolved  than  in  the  previous 
figures,  corresponding  to  a  As  .  °f  0.0030.  One 

striking  difference  between  these  two  cases  is 
near  the  trailing  edge  wake  area.  The  adaptive 
boundary  derivative  matching  term  f^  is  included 

in  the  weighting  function  F(s)  in  Figures  8a  and 
9a,  but  is  not  included  in  Figures  8b  and  9b. 

With  f^  turned  on  (Figure  8b),  the  large  density 

solution  errors  (and  presumably  all  other  solution 
errors)  present  in  the  trailing  edge  zone  of  Fig¬ 
ure  8a  appear  to  reduce  to  the  level  of  the  origi¬ 
nal  fixed  grid  solution  (Figure  3a).  A  disadvan¬ 
tage  of  function  can  be  seen  in  the  adapted 

grid  of  Figure  8b,  however.  The  addition  of 

into  F(s)  seems  to  shear  the  cells  upstream  of  the 
trailing  edge  (i.e.,  near  the  boundaries  of  the 
adaptive  domain),  more  prominently  on  the  lower 
surface  of  the  airfoil.  Perhaps  this  can  be  cor¬ 
rected  by  using  a  modified  form  of  the  function 


In  an  attempt  to  validate  both  the  location 
and  strength  of  the  shocks  predicted  on  the  adapt¬ 
ed  grids  of  Figure  8,  the  resulting  Cp-curves  from 
these  grids  were  compared  with  data  recently  gen¬ 
erated  by  Coakley1*.  Coakley  has  applied  an  im¬ 
plicit  second-order  upwind  scheme  to  the  Euler 
equations  for  identical  flow  conditions  past  the 
NACA0012  airfoil.  The  C-type  grid  used  in  that 
study  was  of  nearly  the  same  dimensions  as  the 
test  grid  of  this  study,  although  grid  point  clus¬ 
tering  near  the  shock  region  was  slightly  higher 
on  Coakley' s  grid. 

Due  to  the  differences  of  the  finite  differ¬ 
ence  structures  employed  in  each  method,  it  was 
anticipated  Chat  an  implicit  upwind  scheme  would 
better  resolve  a  shock  than  an  ADI -type  scheme  on 
similar  grids.  The  curves  in  Figure  10a  show  this 
to  be  the  case.  As  expected,  a  well-defined  shock 
is  observed  in  Coakley's  data,  but  not  in  the  data 
from  the  original  grid.  The  Cp-curves  correspond¬ 
ing  to  the  adaptive  grids  of  Figure.  8,  however, 
more  closely  resemble  the  data  from  the  upwind 

scheme.  In  fact,  in  Figure  10c,  (As  .  =  0.0025), 

min 

the  Cp-curves  match  remarkably  well,  particularly 
in  shock  locations,  which  differ  by  as  little  as 
two  per  cent.  There  is  also  good  Cp  agreement  in 
Figure  10c  downstream  of  the  shock  and  on  the  low¬ 
er  airfoil  surface. 

The  most  disconcerting  region  of  the  adaptive 
grid  Cp  curve  spans  from  the  leading  edge  to  the 
shock  on  the  upper  surface.  The  discrepancy  in  Cp 
values  indicates  that  the  adaptive  grid  Euler  flow 
above  the  airfoil  does  not  accelerate  to  the  Mach 
number  realized  with  the  upwind  scheme.  Assuming 
that  this  problem  was  due  to  truncation  errors  in¬ 
duced  by  insufficient  point  clustering  at  the 
leading  edge,  an  additional  adaptive  grid  run  was 

made,  with  As  .  =  0.0025,  and  with  function  f _ . 

min  2 

Adaption  to  grid  curvature  was  included  to  bring 
more  points  to  the  leading  edge.  As  seen  in  Fig¬ 
ure  lOd,  with  higher  leading  edge  clustering,  the 
Cp  curve  more  closely  matches  Coakley's  data,  al¬ 
though  there  is  an  overshoot  in  pressure  just  be¬ 
fore  the  shock  on  the  newly  adapted  grid.  Note 
also  the  difference  in  the  trailing  edge  Cp  curves 
of  Figure  10c  and  lOd.  This  is  attributable  to 
function  f^,  which  is  used  only  in  Figure  10c. 


Compared  with  the  fixed  grid  solution  of  Figure 
10a,  however,  the  adaptive  grid  solutions  of  Fig¬ 
ures  10b‘10d  more  clearly  resolve  the  upper  sur¬ 
face  midchord  shock  wave. 

RAE2822  Airfoil  Results 

To  further  validate  the  adaption  grid  scheme, 
it  is  desirable  to  compare  numerical  adaptive  grid 
results  with  empirical  data.  Consequently,  numer¬ 
ical  solutions  were  obtained  for  flow  conditions 
equivalent  to  those  presented  in  an  empirical 
study  by  Cook  et  al.lV  In  that  experiment,  ex¬ 
tensive  boundary  layer,  wake  and  pressure  measure¬ 
ments  were  made  for  transonic  flow  past  an  RAE2822 
airfoil.  The  particular  case  selected  for  numeri¬ 
cal  comparison  corresponded  to  turbulent  steady 
flow,  with  a  Mach  number  M  =.73,  an  angle  of  at¬ 
tack  a=3.19°  and  a  Reynolds  number  equal  to  6.52 
million.  Experimentally,  the  boundary  layer  was 
tripped  at  a  distance  of  0.03  chord  lengths  down¬ 
stream  of  the  leading  edge.  Under  these  condi¬ 
tions,  the  boundary  layer  did  not  separate  from 
the  airfoil,  and  a  weak  midchord  shock  formed  on 
the  airfoil  upper  surface.  Considering  these 
facts  and  the  high  Reynolds  number  of  the  flow,  it 
was  reasonable  to  assume  that  the  flow  could  be 
modeled  by  an  inviscid  approximation. 

For  this  reason,  numerical  solutions  were 
again  generated  from  Tassa’s  Navier-Stokes  code  in 
the  Euler  mode.  Except  for  the  angle  of  attack, 
which  was  reduced  to  a=2.57°  from  a  wall  interfer¬ 
ence  correction  formula  suggested  by  Cook  et  al., 
and  except  for  the  Reynolds  number,  flow  condi¬ 
tions  used  were  identical  to  those  jf  the  experi¬ 
ment.  The  inner  detail  of  the  initial  grid  used 
for  these  results  is  presented  in  * igure  11a. 

Only  the  steady-state  solution  was  of  interest,  so 
the  solution  was  advanced  from  impulsive  free 
stream  conditions  using  variable  time  steps,  and 
was  seen  to  converge  after  several  hundred  inte¬ 
grations.  The  density  contours  for  this  case, 
presented  in  Figure  12a,  give  no  indication  of  any 
shock  formation.  With  hopes  of  defining  a  shock, 
the  flow  field  was  solved  again  l..om  an  impulsive 
start,  this  time  with  the  adaptive  grid  solver  em¬ 
ployed  after  every  20  iterations.  The  minimum 

spacing  constraint  was  set  at  As  .  =0.005,  and 

min 

both  functions  and  were  turned  off.  The 

converged  adaptive  grid  for  this  case  is  shown  in 

Figure  lib.  The  large  As  .  constraint  selected 
min 

here  prohibits  grid  point  clustering  to  the  degree 
seen  in  previous  NACA0012  cases.  Nevertheless, 
aside  from  the  leading  edge  region,  the  highest 
point  density  appears  to  be  just  downstream  of  the 
upper  surface  raidchord  region.  As  before,  this 
high  clustering  region  corresponds  to  a  shock 
wave,  seen  in  Figure  12b. 

The  numerical  and  experimental  results  of  this 
flow  are  presented  together  in  Figure  13.  The 
first  of  these  figures  compares  empirical  surface 
pressures  with  numerical  surface  pressures  ob¬ 
tained  on  the  original  grid.  Outside  the  first  20 
per  cont  of  the  airfoil,  the  empirical  and  numeri¬ 
cal  curves  do  not  compare  very  well.  Figure  13b 
compares  the  experimental  data  with  the  solution 
obtained  on  the  adapted  grid  of  Figure  12.  The  Cp 
curves  compare  very  weil  along  the  entire  lower 
surface  of  the  airfoil  In  addition,  the  loca¬ 
tions  of  the  shock  are  in  very  near  agreement.  As 
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is  often  observed,  the  shock  location  in  the  Euler 
solution  is  downstream  of  the  experimental  loca¬ 
tion,  although  the  distance  here  is  rather  small, 
less  than  10  per  cent  of  one  chord  length. 

Once  again,  the  largest  discrepancy  in  the  Cp 
curves  is  on  the  airfoil  upper  surface  upstream  of 
the  shock.  No  attempt  was  made  to  run  another 
case  with  induced  clustering  at  the  leading  edge 
as  was  done  for  the  NACA0012  airfoil  in  the  last 
section.  One  additional  source  of  error  may  have 
been  in  the  effective  angle  of  attack  used  compu¬ 
tationally.  It  is  conceivable  that  better  results 
would  be  realized  if  the  lift  coefficients  C 

L 

rather  than  the  angles  of  attack  are  matched  in 
the  numerical  and  experimental  comparisons.  De¬ 
spite  this  one  region  of  poor  Cp  correlation,  the 
adaptive  grid  routine  has  indeed  clustered  points 
sufficiently  to  resolve  the  shock  wave  seen  in  the 
exper i mental  data  but  not  seen  with  the  fixed 
grid. 

DISCUSSION 

The  practicality  of  a  1 -dimens iona 1  adaptive 
grid  scheme  applied  to  a  2-di mens  tonal  transonic 
airfoil  problem  has  been  demonstrated  in  this 
work.  Advantages  of  a  1 -D  scheme  are  numerous. 
Besides  being  easy  to  formulate  and  easy  to  under¬ 
stand,  this  algorithm  can  be  attached  to  an  exis¬ 
tent  code  with  relatively  few  problems.  Addition¬ 
ally.  the  added  computational  effort  needed  for 
the  scheme  is  minor.  One  run  through  the.  adaption 
subroutine  took  loss  CPU  time  than  one  half  of  one 
time  integration  step.  With  the  adaptive  routine 
utilized  after  every  20  to  30  integration  steps, 
the  increase  in  CPU  time  was  only  about,  two  per 
cent.  Furthermore,  solution  convergence  rates 
with  and  without  adaption  we.re  nearly  the  same 
when  a  solution  was  run  with  a  fixed  At  time  step. 

Unfortunately,  there  also  disadvantages  of  a 
1 -D  adaption  scheme  on  a  2-D  problem.  The  fore¬ 
most  problem  is  that  in  the  equid istr ibution  for¬ 
mulation,  the  final  grid  point  spacing  along  one  r\ 
=  constant  curve  is  almost  totally  independent  of 
the  adaptive  distribution  of  every  other  curve. 

As  a  result,  grid  intersections  are  rarely  orthog¬ 
onal,  and  sometimes  become  so  skewed  that  the  so¬ 
lution  diverges.  Grid  skewness  is  evident  in  Fig¬ 
ure.  3b,  and  to  some  extent ,  in  Figure  8a.  Near 
the  airfoil  leading  edge,  where  arc  curvature  is 
at  its  highest,  skewness  often  became  a  problem. 
With  the  skewness  came  large  solution  errors, 
which  induced  even  more  skewness,  due  to  the  form 
of  the  selected  weighting  function  F(s).  It  was 
almost  always  necessary  to  use  conservative 

(large)  values  of  As  to  insure  a  convergent 
mm  ° 

adaptive  grid.  Because  of  this  inherent  skewness 
problem,  a  technique  which  enforces  orthogona 1 ity 
at  the  grid  intersection  would  bo  highly  desira¬ 
ble  . 

Finally,  as  mentioned  earlier,  grid  speeds  in 
the  transformed  equations  of  motion  were  set  equal 
to  zero  in  this  study,  and  as  a  result,  the  solu¬ 
tion  vector  was  interpolated  after  each  adaptive 
sweep.  Neglecting  grid  speeds  prevented  time-ac¬ 
curate  solutions  from  being  obtained,  however.  If 
the  grid  speeds  were  indeed  calculated  from  back¬ 
ward  differences  and  included  in  the  equations  of 
motion,  time  accurate  convergence  rate  studies 
could  be  made,  and  the  true  effect  of  the  adaptive 
scheme  on  total  computational  time  could  be 
determined . 
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Fig,  2  Test  grid  (  -0.2<x<1.2  ) 


Fig.  3a  Density  contours  of  converged  solution 
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Fig«  3b  Teat  grid  adapted  to  density  gradient  (£^ ) 
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Fig.  10  NACA0012  airfoil  Cp-distribut iom 
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Fig.  11  RAE2822  airfoil  grids 


Original  grid 


Fig,  12  RAE2822  airfoil  density  contours 


Fig,  13  RAE2822  airfoil  Cp-distributions 


